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Algebra of Vectors and Tensors 

Whereas heat and mass are scalars, fluid mechanics concerns transport of momentum, which is a 
vector. Heat and mass fluxes are vectors, momentum flux is a tensor. Consequently, the mathematical 
description of fluid flow tends to be more abstract and subtle than for heat and mass transfer. In an 
effort to make the student more comfortable with the mathematics, we will start with a review of the 
algebra of vectors and an introduction to tensors and dyads. A brief review of vector addition and 
multiplication can be found in Greenberg,* pages 132-139. 

Scalar - a quantity having magnitude but no direction (e.g. temperature, density) 

Vector - (a.k.a. 1st rank tensor) a quantity having magnitude and direction (e.g. velocity, force, 
momentum) 

(2nd rank) Tensor - a quantity having magnitude and two directions (e.g. momentum flux, 
stress) 


a 


Vector Multiplication 

Given two arbitrary vectors a and b, there are three types of vector products 
are defined: 


Notation Result Definition 


Dot Product 

ab 

scalar ab cos0 

Cross Product 

axb 

vector ab | sin0 | n 


where 0 is an interior angle (0 < 0 < k) and n is a unit vector which is normal to both a and b. 
sense of n is determined from the "right-hand-rule"'* 

Dyadic Product ab tensor 



b 


The 


* Greenberg, M.D., Foundations Of Applied Mathematics, Prentice-Hall, 1978. 

♦ The “ right-hand rule”: with the fingers of the right hand initially pointing in the direction of the first 
vector, rotate the fingers to point in the direction of the second vector; the thumb then points in the 
direction with the correct sense. Of course, the thumb should have been normal to the plane containing 
both vectors during the rotation. In the figure above showing a and b, axb is a vector pointing into the 
page, while bxa points out of the page. 


Copyright © 2000 by Dennis C. Prieve 



06-703 


2 


Fall, 2000 


In the above definitions, we denote the magnitude (or length) of vector a by the scalar a. Boldface will 
be used to denote vectors and italics will be used to denote scalars. Second-rank tensors will be 
denoted with double-underlined boldface; e.g. tensor J. 


Definition of Dyadic Product 

Reference: Appendix B from Happel & Brenner. v The word “dyad” comes from Greek: “dy” 
means two while “ad” means adjacent. Thus the name dyad refers to the way in which this product is 
denoted: the two vectors are written adjacent to one another with no space or other operator in 
between. 

There is no geometrical picture that I can draw which will explain what a dyadic product is. It's best 
to think of the dyadic product as a purely mathematical abstraction having some very useful properties: 

Dyadic Product ab - that mathematical entity which satisfies the following properties (where a, 
b, v, and w are any four vectors): 

1. ab • v = a(b • v) [which has the direction of a; note that ba • v = b(a • v) which has the direction of 
b. Thus ab ± ba since they don’t produce the same result on post-dotting with v.] 

2. v • ab = (v • a)b [thus v • ab ^ ab • v] 

3. abxv = a(bxv) which is another dyad 

4. vxab = (vxa)b 

5. ab : vw = (a • w)(b • v) which is sometimes known as the inner-outer product or the double-dot 
product .* 

6. a(v+w) = av+aw (distributive for addition) 

7. (v+w)a = va+wa 

8. (.v+/)ab = sab+/ab (distributive for scalar multiplication-also distributive for dot and cross 
product) 

9. sab = (,sa)b = a(.vb) 


v Happel, J., & H. Brenner, Low Reynolds Number Hydrodynamics, Noordhoff, 1973. 

* Brenner defines this as (a • v)(b • w). Although the two definitions are not equivalent, either can be 
used - as long as you are consistent. In these notes, we will adopt the definition above and ignore 
Brenner's definition. 
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Decomposition into Scalar Components 

Three vectors (say e ] , e2, and 63) are said to be linearly independent if none can be expressed 
as a linear combination of the other two (e.g. i, j, and k). Given such a set of three LI vectors, any 
vector (belonging to E^) can be expressed as a linear combination of this basis : 


V — + VoC 2 P3C3 


where the Vj are called the scalar components of v. Usually, for convenience, we choose 
orthonormal vectors as the basis: 


e L e / = 5 // = 



although this is not necessary. 5 ;/ - is called the Kronecker delta. Just as the familiar dot and cross 
products can written in terms of the scalar components, so can the dyadic product: 


vw — (v ^ +v 26 2+v 363X1+ 36 j+hd 6 2+1+36 3) 


= (vie 1 Xw 1 e 1 )+(v 1 e 1 )(w 2 e2)+ ... 


= v 1 w 1 e 1 e 1 +v 1 i+2e 1 e2+ ... (nine terms) 

where the e,-ey are nine distinct unit dyads. We have applied the definition of dyadic product to 
perform these two steps: in particular items 6, 7 and 9 in the list above. 

More generally any nth rank tensor (in E^) can be expressed as a linear combination of the J l unit 11- 
ads. For example, if n= 2, 3 /7 =9 and an n-ad is a dyad. Thus a general second-rank tensor can be 
decomposed as a linear combination of the 9 unit dyads: 

1 = 7 ’ll e l e l +7 ’l2 e l e 2+ ••• = ^i=l,3^j=l,3 T ifi e j 

Although a dyad (e.g. vw) is an example of a second-rank tensor, not all 
2nd rank tensors T can be expressed as a dyadic product of two vectors. 

To see why, note that a general second-rank tensor has nine scalar 
components which need not be related to one another in any way. By 
contrast, the 9 scalar components of dyadic product above involve only six 
distinct scalars (the 3 components of v plus the 3 components of w). 

After a while you get tired of writing the summation signs and lim its. So an 
abbreviation was adopted whereby repeated appearance of an index implies summation over the three 
allowable values of that index: 



Tee 

U'J 
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This is sometimes called the Cartesian (implied) summation convention. 


Scalar Fields 

Suppose I have some scalar function of position (. x,y,z ) which is continuously differentiable , that 
is 


f=f(x,y,z) 

and df/dx, df/dy, and df/dz exist and are continuous throughout some 3-D region in space. This 
function is called a scalar field. Now consider / at a second point which is differentially close to the 
first. The difference in / between these two points is 
called the total differential of/: 

f(x + dx, y + dy, z.+dz) -f(x,y,z) = df 

For any continuous function f(x,y,z), df is linearly related 
to the position displacements, dx, dy and dz. That 
linear relation is given by the Chain Rule of 
differentiation: 

df df df 

df = — dx -\ dy dz 

dx dy dz 

Instead of defining position using a particular coordinate 
system, we could also define position using a position vector r: 

r = .ri + yj + zk 

The scalar field can be expressed as a function of a vector argument, representing position, instead of a 
set of three scalars: 



f=m 

Consider an arbitrary displacement away from the point r, which we denote as <r/r to emphasize that the 
magnitude | dr | of this displacement is sufficiently small that /( r) can be linearized as a function of 
position around r. Then the total differential can be written as 
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df=f(r + dr)-f(r) 


Gradient of a Scalar 

We are now is a position to define an important vector associated 
with this scalar field. The gradient (denoted as Vf) is defined 
such that the dot product of it and a differential displacement 
vector gives the total differential: 


/ 



M 


df = dr -Vf 


EXAMPLE: Obtain an explicit formula for calculating the gradient in Cartesian* coordinates. 

Solution : r = xi + yj + zk 

r+dr = (x+d.x)i + (y+dy)} + (z+dz)k 
subtracting: dr = (dx)i + (c/y)j + (dz) k 

V/ = (V/) r i + (V/) y j + (V/) 7 k 
dr- Vf = [(dx)i + ...]• [(Vf) x i + ...] 

df = ( V f) x dx + (Vf) y dy + ( Vf)Jz ( 1 ) 

Using the Chain rule: df = (df/r)x)dx + ( df/dy)dy + ( df/dz)dz (2) 

According to the definition of the gradient, (1) and (2) are identical. Equating them and collecting terms: 

l(Vf) x -(df/dx)\dx + [(Vf) y -(df/dy)]dy + [(V f) z -(df/dz)]dz = 0 

Think of dx, dy, and dz as three independent variables which can assume an i nfi nite number of values, 
even though they must remain small. The equality above must hold for all values of dx, dy, and dz. The 
only way this can be true is if each individual term separately vanishes:** 


*Named after French philosopher and mathematician Rene Descartes (1596-1650), pronounced "day- 
cart", who first suggested plotting fix) on rectangular coordinates 

** For any particular choice of dx, dy, and dz, we might obtain zero by cancellation of positive and 
negative terms. However a small change in one of the three without changing the other two would cause 
the sum to be nonzero. To ensure a zero-sum for all choices, we must make each term vanish 
independently. 
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So (Vf) x = df/dx, (V/) y = df/dy, and (V// = df/dz, 

i • v df.df.df 

leavmg V/ = — 1 + — j + — k 

dx dy c)z 

Other ways to denote the gradient include: 

V/= grad/= df/dr 

Geometric Meaning of the Gradient 

1) direction: V/( r) is normal to the f=c onst surface passing through the point r in the direction of 
increasing/. V/ also points in the direction of steepest ascent of/. 

2) magnitude: |V/ is the rate of change of / with 
distance along this direction 

What do we mean by an '/= const surface"? Consider an 
example. 

Example : Suppose the steady state temperature profile 
in some heat conduction problem is given by: 

T(x,y,z) = x 2 + y 2 + z 2 

Perhaps we are interested in VT at the point (3,3,3) 
where 7=27. VT is normal to the 7=const surface: 

x 2 + y 2 + z 2 = 27 

which is a sphere of radius -Jr] .* 

Proof of 1 ). Let's use the definition to show that these geometric meanings are correct. 

df= dv Vf 


* A vertical bar in the left margin denotes material which (in the interest of time) will be omitted from the 
lecture. 
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Consider an arbitrary /. A portion of the /=const surface 
containing the point r is shown in the figure at right. Choose a 
dr which lies entirely on/=const. In other words, the surface 
contains both r and r+dr, so 

/( r) =/( r+dr) 

and df = /(r+dr)-/(r) = 0 

Substituting this into the definition of gradient: 

df= 0 = dr- V/ = Ur | I V/' | cosG 



Since I dr | and I V/ 1 are in general not zero, we are forced 

to the conclusion that cos0=O or 0=90° . This means that V/ is normal to dr which lies in the surface. 


2) can be proved in a similar manner: choose dr to be parallel to V/ Does V/ point toward higher or 
lower values of/? 


Applications of Gradient 

• find a vector pointing in the direction of steepest ascent of some scalar field 

• determine a normal to some surface (needed to apply b.c.’s like n - v = 0 for a boundary which is 
impermeable) 

• determine the rate of change along some arbitrary direction: if n is a unit vector pointing along some 
path, then 


- v/ =! 


is the rate of change of / with distance (s) along this path given by n. df /ds is called the directed 
derivative of/. 


Curvilinear coordinates 

In principle, all problems in fluid mechanics and transport could be solved using Cartesian 
coordinates. Often, however, we can take advantage of symmetry in a problem by using another 
coordinate system. This advantage takes the form of a reduction in the number of independent variables 
(e.g. PDE becomes ODE). A fa mil iar example of a non-Cartesian coordinate system is: 
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Cylindrical Coordinates 

r = (. x 2 +y 2 ) 1/2 v = rc osO 

0 = tan ~ l (y/x) y = rsinO 

z = z z = z 

Vectors are decomposed differently. Instead of 
inR.C.C.S.: v = v^i + v-yj + v z k 

in cylindrical coordinates, we write 
in cyl. coords.: v = v r e r + vqCq + v z e z 

where e r , eg, and e, are new unit vectors pointing the r, 0 and z directions. We also have a different 
set of nine unit dyads for decomposing tensors: 

e r e r , e r e 0 , e r e z , eee r , etc. 

Lik e the Cartesian unit vectors, the unit vectors in cylindrical coordinates form an orthonormal set of 
basis vectors for E 3 . U nlik e Cartesian unit vectors, the orientation of e r and eg depend on position. In 
other words: 

e r = e r ( 0 ) 
e 9 = e e ( 0 ) 
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Spherical Coordinates 


z 


0 . 


0 . 


0 . 


0 , 



Spherical coordinates ) are defined relative to Cartesian coordinates as suggested in the 
figures above (two views of the same thing). The green surface is the xy-plane, the red surface is the 
xz-plane, while the blue surface (at least in the left image) is the yz-plane. These three planes intersect at 
the origin (0,0,0), which lies deeper into the page than (1,1,0). The straight red line, drawn from the 
origin to the point (r.G.p) 4- has length r, The angle 0 is the angle the red line makes with the z-axis (the 
red circular arc labelled 0 has radius r and is subtended by the angle 0). The angle (f> (measured in the 
xy-plane) is the angle the second blue plane (actually it’s one quadrant of a disk) makes with the xy- 
plane (red). This plane which is a quadrant of a disk is a (|)=const surface: all points on this plane have 
the same (|> coordinate. The second red (circular) arc labelled <|) is also subtended by the angle <|) . 


* This particular figure was drawn using r = 1,0 = tc/ 4 and (J> = tc/3. 
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A number of other (|)=const planes are 
shown in the figure at right, along with a 
sphere of radius r=l . All these planes 
intersect along the z-axis, which also passes 
through the center of the sphere. 


a: = r sinOcos (|) 
y = rsinOsint]) 

z = rcos 0 


The position vector in spherical coordinates 
is given by 

r = vi+yj+zk = r e r (0,<))) 



0 = tan 


(]) = tan (y/.r) 


In this case all three unit vectors depend on 
position: 


e r = e r (0,<>), e e = e 0 (9,<|)), and e^= e^) 


where e r is the unit vector pointing the direction of increasing r, holding 0 and (|) fixed; e e is the unit 
vector pointing the direction of increasing 0, holding r and <|> fixed; and e 0 is the unit vector pointing the 
direction of increasing (|> , holding r and 0 fixed. 


These unit vectors are shown in the figure at right. 
Notice that the surface (|)=const is a plane containing the 
point r itself, the projection of the point onto the xy-plane 
and the origin. The unit vectors e r and e e lie in this plane 


as well as the Cartesian 



unit circle on 
<(> = const 
surface 


unit vector k 
denoted e,). 


(sometimes 



If we tilt this (|)=const plane 

into the plane of the page (as in the sketch at left), we can more easily see 
the relationship between these three unit vectors: 

e~ =(cos0)e,. -(sinO)eQ 

This is obtained by determined from the geometry of the right triangle in 
the figure at left. When any of the unit vectors is position dependent, we 
say the coordinates are: 


Copyright © 2000 by Dennis C. Prieve 


06-703 


11 


Fall, 2000 


curvilinear - at least one of the basis vectors is position dependent 

This will have some profound consequences which we will get to shortly. But first, we need to take 
“time-ouf ’ to define: 


Differentiation of Vectors w.r.t. Scalars 

Suppose we have a vector v which depends on the scalar parameter t : 

v = v(t) 

For example, the velocity of a satellite depends on time. What do we mean by the “derivative” of a 
vector with respect to a scalar. As in the Fundamental Theorem of Calculus, we define the derivative 
as: 


d\ _ lim J \(t + At) - v(t) 
dt Ar— >o\ At 


Note that d\/dt is also a vector. 

EXAMPLE: Compute deddQ in cylindrical coordinates. 


Solution : From the definition of the derivative: 

de r _ lim je,.(0 + A0) 
dQ A0— >0 [ A0 


Since the location of the tail of a vector is not part 
of the definition of a vector, let's move both 
vectors to the origin (keeping the orientation 
fixed). Using the parallelogram law, we obtain the 
difference vector. Its magnitude is: 

|e r (0+A0)-e r (0)| = 2sin4f- 

Its direction is parallel to ee(0+A0/2), so: 

e r (0 + A0) -e r (0) = 

Recalling that si nr tends to x as x— >0, we have 
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lim {e r (6 + A0) -e r (0)} = A0e e (0) 
A0— >0 

Dividing this by A0, we obtain the derivative: 

r/e/r/0 = eg 

Similarly, de^/dQ = -e r 

One important application of “differentiation with respect to a 
scalar” is the calculation of velocity, given position as a function of 
time. In general, if the position vector is known, then the velocity 
can be calculated as the rate of change in position: 

r=r (0 

v = dr/dt 

Similarly, the acceleration vector a can be calculated as the 
derivative of the velocity vector v: 

a = d\/dt 

EXAMPLE: Given the trajectory of an object in 
cylindrical coordinates 

r = r(t), 0 = 0(0, and z = z(t ) 

Find the velocity of the object. 

Solution : First, we need to express r in in terms of the 
unit vectors in cylindrical coordinates. Using the figure at 
right, we note by inspection that* x 

r (r,0,z) = re r (0) + ze. 

Now we can apply the Chain Rule: 


*Recalling that r = xi + _vj + zk in Cartesian coordinates, you might be tempted to write r = re r + 0e e + 
ze z in cylindrical coordinates. Of course, this temptation gives the wrong result (in particular, the units 
of length in the second term are missing). 
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Jr = 


'dr ' 


Jr + 


'Or 

7e r (9) 

1 #3 
e e 


JO + 


Jr 

ife 


dz = e r Jr + reg J0 + e z dz 


Dividing by dt, we obtain the velocity: 


_ Jr _ Jr(?) J0(t) Jz(f) 

0 ^ I /* Gq I 0 -- 

dt 1^3 1 2 l B 1^3 

*7 v e v z 


Vector Fields 

A vector field is defined just like a scalar field, except that it's a vector. Namely, a vector field is a 
position-dependent vector: 


v = v(r) 

Common examples of vector fields include force fields, like the gravitational force or an electrostatic 
force field. Of course, in this course, the vector field of greatest interest is: 


Fluid Velocity as a Vector Field 

Consider steady flow around a submerged object. What do we mean by “fluid velocity?” There 
are two ways to measure fluid velocity. First, we could add tracer particles to the flow and measure the 
position of the tracer particles as a function of time; differentiating position with respect to time, we 
would obtain the velocity. ♦ A mathematical “tracer particle” is called a “material point:” 

Material point - fluid element - a given set of fluid molecules whose location may change with 
time.* 


♦ Actually, this only works for steady flows. In unsteady flows, pathlines, strea klin es and streamlines 
differ (see “Streamlines, Pathlines and Streaklines” on page 65). 

* In a molecular-level description of gases or liquids, even nearby molecules have widely different 
velocities which fluctuate with time as the molecules undergo co lli sions. We will reconcile the 
molecular-level description with the more common continuum description in Chapter 4. For now, we 
just state that by “location of a material point” we mean the location of the center of mass of the 
molecules. The “point” needs to contain a statistically large number of molecules so that if/) converges 
to a smooth continuous function. 
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Suppose the trajectory of a material point is given by: 


r=r (t) 
clr 

Then the fluid velocity at any time is v = — 


(3) 


A second way to measure fluid velocity is similar to the “bucket-and-stopwatch method.” We measure 
the volume of fluid crossing a surface per unit time: 


n- v = lim 

A Cl — >0 

where A a is the area of a surface element having a unit 
normal n and Ac/ is the volumetric flowrate of fluid crossing 
A a in the direction of n. 



v 



When A a is small enough so that this quotient has converged 

in a mathematical sense and A a is small enough so that the surface is locally planar so we can denote its 
orientation by a unit normal n, we can replace A a by da and Ac/ by dq and rewrite this definition as: 

dq = n • v da (4) 

This is particularly convenient to compute the 
volumetric flowrate across an arbitrary curved 
surface, given the velocity profile. We just have to 
sum up the contribution from each surface element: 

Q = J n -v da 

A 


Partial & Material Derivatives 



Let f=f(r,t) 

represent some unsteady scalar field (e.g. the unsteady temperature profile inside a moving fluid). There 
are two types of time derivatives of unsteady scalar fields which we will find convenient to define. In the 
example in which/ represents temperature, these two time derivatives correspond to the rate of change 
(denoted genetically as df Idt) measured with a thermometer which either is held stationary in the 
moving fluid or drifts along with the local fluid. 

partial derivative - rate of change at a fixed spatial point: 
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di y dt J (lr=0 

where the subscript <r/r=0 denotes that we are evaluating the derivative along a path* on 
which the spatial point r is held fixed. In other words, there is no displacement in 
position during the time interval dt. As time proceeds, different material points occupy 
the spatial point r. 

material derivative (a.k.a. substantial derivative) - rate of change within a particular material 
point (whose spatial coordinates vary with time): 

£l = (dC 

Dt ^ dt J c i r=V( j t 

where the subscript dr = v dt denotes that a displacement in position (corresponding to 
the motion of the velocity) occurs: here v denotes the local fluid velocity. As time 
proceeds, the moving material occupies different spatial points, so r is not fixed. In 
other words, we are following along with the fluid as we measure the rate of change of 

A relation between these two derivatives can be derived using a generalized vectorial form of the Chain 
Rule. First recall that for steady (independent of t) scalar fields, the Chain Rule gives the total 
differential (in invariant form) as 

df - f{ r + dr)- f(r) = dr- Vf 

When l is a variable, we just add another contribution to the total differential which arises from changes 
in t, namely dt. The Chain Rule becomes 

df = f(r + dr,t + dt)-f(r,t) =^-dt + dr- V/ 

The first term has the usual Chain-Rule form for changes due to a scalar variable; the second term gives 
changes due to a displacement in vectorial position r. Dividing by dt holding R fixed yields the material 
derivative: 


* By “path” I mean a constraint among the independent variables, which in this case are time and 
position (e.g. x,y,z and 1). For example, I might vary one of the independent variables (e.g. x) while 
holding the others fixed. Alternatively, I might vary one of the independent variables (e.g. t) while 
prescribing some related changes in the others (e.g. x(t), y(t ) and z.(t)). In the latter case, I am 
prescribing (in parametric form) a trajectory through space, hence the name “path.” 
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PI = , 

Dt 


( ,ir \ 


df_ 

dt 


)dr=\dt 


dt 


^ dt \ 




+ 


<r/r 

dt 


•v/ 


But (dr/dl) is just v, leaving: 


d/ a/ 

— = — + v Vf 
Dt dt 


This relationship holds for a tensor of any rank. For example, the material derivative of the velocity 
vector is the acceleration a of the fluid, and it can be calculated from the velocity profile according to 


Dv dv 

a = = — + v-Vv 

Dt dt 


We will define Vv in the next section. 


Calculus of Vector Fields 

Just like there were three kinds of vector multiplication which can be defined, there are three kin ds 
of differentiation with respect to position. 


Shortly, we will provide explicit definitions of these 


Notation 

Result 

quantities in terms of surface integrals. Let me 
introduce this type of definition using a more familiar 

Divergence 

V- v 

scalar 

quantity: 

Curl 

Vxv 

vector 


Gradient 

Vv 

tensor 


Gradient of a Scalar (Explicit) 


Recall the previous definition for gradient: 

/ =/ (r): df = dr • Vf (implicit def n of Vf ) 

Such an implicit definition is like defining f' (x) as that function associated with fix) which yields: 

f=f(x): df = (dx) f ' (implicit def n off ’ ) 

An equivalent, but explicit, definition of derivative is provided by the Fundamental Theorem of the 
Calculus: 


f' {x \= lim [ /(* + Ax) -/(*) ] _ df 
Ax — > 0 \ Ax j dx 

We can provide an analogous definition of Vf 


(explicit def n off) 
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V/ = 


lim 

V-»0 


V 



(explicit def n of V/ ) 


where/ = any scalar field 

A = a set of points which constitutes any 
closed surface enclosing the point r 
at which V/ is to be evaluated 

V = volume of region enclosed by A 

da = area of a differential element (subset) of 
A 

n = unit normal to da, pointing out of region 
enclosed by A 



lim (V— >0) = lim it as aH dimensions of A shrink to zero (in other words, A collapses about the 
point at which V/ is to be defined.) 


What is meant by this surface integral? Imagine A to be the skin of a potato. To compute the integral: 

1) Carve the skin into a number of elements. Each element must be sufficiently small so that 


• element can be considered planar (i.e. n is practically constant over the element) 

• /is practically constant over the element 


2) For each element of s kin , compute n f da 

3) Sum yields integral 


This same type of definition can be used for each of the three spatial derivatives of a vector field: 


Divergence, Curl, and Gradient 


Divergence 


V- v 


lim 


V -> oj 


V 


Jn -\da > 

A 


Curl 


V x v 


lim 

V — » 0 


J_ 

V 


J n x \da > 

A 
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Gradient 


Vv = 


lim 

V^O 


V 


J ny da > 

A 


Physical Interpretation of Divergence 


Let the vector field v = v(r) represent the steady-state velocity profile in some 3-D region of space. 
What is the physical meaning of V • v? 


• n • v da = dq = volumetric flowrate out through da (cm 3 /s). This quantity is positive for 
outflow and negative for inflow. 

• l A n • v da = net volumetric flowrate out of enclosed volume (crnVs). This is also positive for 
a net outflow and negative for a net inflow. 

• (1/V)Ja n • v da = flowrate out per unit volume (s' 1 ) 


• V • V 


> 0 for an expanding gas (perhaps T T or pi) 

< = 0 for an incompressible fluid (no room for accumulation) 
< 0 for a gas being compressed 


• V • v = volumetric rate of expansion of a differential element of fluid per unit volume of that 
element (s' 1 ) 


Calculation of V • v in R.C.C.S. 

Given: v = v x (x,y,z) i + v y (x,y,z ) j + vfx,y,z) k 

Evaluate V • v at (. x 0 ,y 0 ,z 0 ). 


Solution : Choose A to be surface of rectangular 
parallelopiped of dimensions Ax,Ay,Az with one comer 

at x 0 ,y 0 ,z 0 . 



So we partition A into the six faces of the parallelopiped. 
The integral will be computed separately over each face: 

|n-v<r/a= |n-v<r/<3+ J n- v da + L + Jn -v da 

A Ay A 2 A 6 

Surface A 1 is the x=x Q face: n = - i 



n - v = -i • v = -v x (x 0 ,y,z) 
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z 0 + Az y a +A y 

J n -\da= | | -v x (x 0 ,y,z)dydz 


Using the Mean Value Theorem: 

where 

and 

Surface A 2 is the x=x a +Ax face: 


y a 

= -v x (x 0 ,y',z)AyAz 

y a ^y'^ y 0 +A y 


Z a < Z < Z 0 +Az 


n= +i 


n • v = i • v = v x (x a +Ax ,y,z) 


z 0 +Az y a +Ay 

| n-v<r/a= J J v X (x 0 + Ax,y,z) dydz 

A \ z o y a 

Using the Mean Value Theorem: = v x (x 0 +Ax,y",z")AyAz 

where y Q < y" < y () +Ay 

and z 0 < z < z 0 +Az 

Hie sum of these two integrals is: 

J +J =[v x (x 0 + Ax,y",z")~v x (x 0 ,y',z')]AyAz 
A\ A 2 


Dividing by V = AxAyAz: 


— f n- v da = 
V . J . 


A l + A 2 


v x (Xq +A X,y",z")-v x (x 0 ,y',z') 
Ax 


Letting Ay and Az tend to zero: 


lim 

Ay,Az — > 0 


— f n- vda 

V . J . 


A l + a 2 


r , ( -V, , + Av, , Zq ) - i’ A ( Xp , y ( , , z„ ) 
Ax 


Finally, we take the lim it as A,r tends to zero: 
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lim 

1 r 

dv v 

< 

— n- vda 

> = x 

V->0 

V J 

dx 


l A l +A 2 J 


Similarly, from the two y=const surfaces, we obtain: 


lim 

1 r 

ro 

II 

A, 

< 

— n- vda 

V-»0 

V J 

A 3 +A 4 

dy | 


y 0 ’Z 


O 


y Q ^o 


and from the two z=const surfaces: 


lim 

1 r 

dv. 

< 

— n- vda 


V-a0 

V J 

dz 


a 5 + a 6 


Summing these three contributions yields the divergence: 


dv x 

V- v = x 


+ ■ 


dvy dv. 


-+• 


dx dy dz 


Evaluation of Vxv and Vv in R.C.C.S. 

In the same way, we could use the definition to determine expressions for the curl and the gradient. 


Vxv 


^dv, dv v \ fdv x dv z V ( dv v dv ^ 


d y dz 


i + 


d z dx 


J + 


dx dy 


Hie formula for curl in R.C.C.S. turns out to be expressible as a determinant of a matrix: 


j k 

d d 


dx dy dz 


v x v y v z 


dv, dv y \ 
dy dz 


1 + 


dv, 


dv^ 


dz dx 




dv v dv , 


dx dy 


But remember that the determinant is just a mnemonic device, not the definition of curl. The gradient 
of the vector v is 


Vv 


I I 

i= 1 7=1 


dv j 

d^T e/Cj ' 
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where x i = x, x 2 = y, and .*3 = z, v \ = v x , etc. 


Evaluation ofV-\, Vxv and Vv in Curvilinear Coordinates 
Ref: Greenberg, pl75 

These surface-integral definitions can be applied to any coordinate system. On HWK #2, we 
obtain V • v in cylindrical coordinates. 

More generally, we can express divergence, curl and gradient in terms of the metric coefficients 
for the coordinate systems. If u,v,w are the three scalar coordinate variables for the curv ili near 
coordinate system, and 

v = x(u,v,w) y = y(u,v,w ) z = z(u,v,w ) 

can be determined, then the three metric coefficients — h\, I 12 and /13 — are given by 

h\(u,v,w) = -Jx^ + yl +zl 

h 2 (u,v,w ) = -Jxy + y\, + 4 

h 3 (u,v,w) = 7*5, + yl> + i 

where letter subscripts denote partial differentials while numerical subscripts denote component, and the 
general expressions for evaluating divergence, curl and gradient are given by 

. . . . Id/ Id/ Id/ 

gradient of scalar: V/ = — — H — eo H 3-^3 

h\ ou hi dv /13 dii’ 


divergence of vector: 


V-v 


1 

h\h 2 h 3 


_d_ 

d u 


{ h 2 h 3 v l ) + 


_d_ 

dv 


i h l h 3 v 2 ) + 


_d_ 

d w 


{ h l h 2 v 3 ) 


curl of vector: 


V x v = 


1 

h\h 2 h 3 



h v l 



h 2 v 2 


^ 3 e 3 



^ 3 V 3 


These formulas have been evaluated for a number of common coordinate systems, including R.C.C.S., 
cylindrical and spherical coordinates. The results are tabulated in Appendix A of BSL (see pages 738- 
741). These pages are also available online: 
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rectangular coords . 
cylindrical coords : 
spherical coords : 


Physical Interpretation of Curl 

To obtain a physical interpretation of Vxv, let’s consider a particularly simple flow field which is 
called solid-body rotation. Solid-body rotation is simply the velocity field a solid would experience if 
it was rotating about some axis. This is also the velocity field eventually found in viscous fluids 
undergoing steady rotation. 


Imagine that we take a container of fluid 
(like a can of soda pop) and we rotate the 
can about its axis. After a transient period 
whose duration depends on the dimensions 
of the container, the steady- state velocity 
profile becomes solid-body rotation. 

A material point imbedded in a solid would 
move in a circular orbit at a constant angular 
speed equal to <T> radians per second. The 
corresponding velocity is most easily 
described using cylindrical coordinates with 
the z-axis oriented perpendicular to the 
plane of the orbit and passing through the 
center of the orbit. Then the orbit lies in a 
z=const plane. The radius of the orbit is the 
radial coordinate r which is also constant. 
Only the 9-coordinate changes with time and 
it increases linearly so that dQ/dt = const = 

a 

In parametric form in cylindrical coordinates, 
the trajectory of a material point is given by 



e, (9(0)) 


r{t) = const, z(t) = const, 0(f) = Qt 


The velocity can be computed using the formulas developed in the example on page 12: 

dQ(t) dz(t 
r+r —eg H — 

123 1 2*3 12(3 

0 rQ. 0 


■e z = rOeg 
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Alternatively, we could deduce v from the definition of derivative of a 
vector with respect to a scalar: 


Dr Ar rd0e 0 d 0 _ 

v = — = lim — = = r — e 0 = ri 2e 0 

Dt a?->oA/ dt dt 


More generally, in invariant form (i.e. in any coordinate system) the 
velocity profile corresponding to solid-body rotation is given by 




x r, 



where Q is called the angular velocity vector and iy, is the position vector* whose origin lies 
somewhere along the axis of rotation. The magnitude of Q is the rotation speed in radians per unit time. 
It’s direction is the axis of rotation and the sense is given by the “right-hand rule.” In cylindrical 
coordinates, the angular velocity is 

Q. = Qe- 

and the position vector is r„ = re, + ze~ 

Taking the cross product of these two vectors (keeping the order the same as in (5)): 

v( r ) = ^ §,- + z&fe f §- = ^e 0 

e e 0 

To obtain this result we have used the fact that the cross product of any two pai 
(because the sine of the angle between them is zero — recall definition of cross 
product on pi). 

The cross product of two distinct unit vectors in any right-handed coordinate 
system yields a vector parallel to the third unit vector with a sense that can be 
remembered using the figure at right. If the cross product of the two unit vectors 
corresponds to a “clockwise” direction around this circle, the sense is positive; in 
a “counter-clockwise” direction, the sense is negative. In this case, we are 
crossing e - with e r which is clockwise; hence the cross product is +e 0 . 

Now that we have the velocity field, let’s compute the curl. In cylindrical coordinates, the formula for 
the curl is obtained from p739 of BSL: 



The subscript “p” was added here to avoid confusing the cylindrical coordinate r with the magnitude 


of the position vector. Note that in cylindrical coordinates, 


=7777 


* r. 
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1 dv z dv e V , r <K dv z 
r d9 dz J \ dz dr 



1 d(n ; e) 

r dr 



Substituting v r = 0 v e = rO v z = 0 

we obtain V x v = 20e : = 20. 

Thus the curl turns out to be twice the angular velocity of the fluid elements. While we have only shown 
this for a particular flow field, the results turns out to be quite general: 

V x v = 20 


Vector Field Theory 

There are three very powerful theorems which constitute “vector field theory:” 

• Divergence Theorem 

• Stokes Theorem 

• Irrotational <=> Conservative ^Derivable from potential 

Divergence Theorem 

This is also known as ‘ Uauss * Divergence Theorem ” or ‘ TJreen’s Formula ” (by Landau & 
Lifshitz). Let v be any (continuously differentiable) vector field and choose A to be any (piecewise 
smooth, orientable) closed surface; then 

j n- v da = J V- v dV 

A V 

where V is the region enclosed by A and n is the 
outward pointing unit normal to the differential 
surface element having area da. Although we will 
not attempt to prove this theorem, we can offer the 

surface A 

* Carl Friedrich Gauss (1777-1855), German mathematician, physicist, 
the greatest mathematician of his time and the equal of Archimedes and Isaac Newton, Gauss made 
many discoveries before age twenty. Geodetic survey work done for the governments of Hanover and 
Denmark from 1821 led him to an interest in space curves and surfaces. 
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following rationalization. Consider the lim it in which all dimensions of the region are very small, i.e. 
V— >0. When the region is sufficiently small, the integrand (which is assumed to vary continuously with 
position)* is just a constant over the region: 

V • v = const, inside V 


r v 

jn-\da = Jv-vr/V =(V-v) J c IV 

A V VV J 


(V-v)V 


Solving for the divergence, we get the definition back (recalling that this was derived for V— >0): 


V 



J n- \da 

A 


Thus the divergence theorem is at least consistent with the definition of divergence. 


Corollaries of the Divergence Theorem 

Although we have written the Divergence Theorem for vectors (tensors of rank 1), it can also be 
applied to tensors of other rank: 


J nf da = J V/ dV 

A V 

jn-xda = J V-xdV 

A V 

One application of the divergence theorem is to simplify the evaluation of surface or volume integrals. 
However, we will use GDT mainly to derive invariant forms of the equations of motion: 

Invariant : independent of coordinate system. 

To illustrate this application, let’s use GDT to derive the continuity equation in invariant form. 


* This is a consequence of v being “continuously differentiable”, which means that all the partial 
derivatives of all the scalar components of v exist and are continuous. 
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The Continuity Equation 

Let p(r,t) and v(r ,t) be the density and fluid velocity. What relationship between them is imposed 
by conservation of mass? 


For any system, conservation of mass means: 

f rate of acc. 1 f net rate of 

[of total mass] [mass entering 

Let's now apply this principle to an arbitrary system 
whose boundaries are fixed spatial points. Note 
that this system, denoted by V can be macroscopic 
(it doesn’t have to be differential). The boundaries 
of the system are the set of fixed spatial points 
denoted as A. Of course, fluid may readily cross 
these mathematical boundaries. 


v 



Subdividing V into many small volume elements: 


dm = pdV 
M = J dm = J p dV 

V 


dM 

dt 


f \ 




where we have switched the order of differentiation and integration. This last equality is only valid if the 
boundaries are independent of 1. Now mass enters through the surface A. Subdividing A into small 
area elements: 


n = outward unit normal 
n • v da = vol. flowrate out through da (cnriVs) 
p(n • v)da = mass flowrate out through da (g/s) 

{ }= f p (n-v) da = f n- ( pv ) da = fv-(pv)<r/V 

[mass leaving J J J J 

The third equality was obtained by applying GDT. Substituting into the general mass balance: 
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-Jv-(pv)c/V 

V 


Since the two volume integrals have the same limits of integration (same domain), we can combine them: 



dV = 0 


Since V is arbitrary, and since this integral must vanish for all V, the integrand must vanish at every 
point:* 


|£+V-(pv) = 0 

which is called the equation of continuity. Note that we were able to derive this result in its most 
general vectorial form, without recourse to any coordinate system and using a finite (not differential) 
control volume. In the special case in which p is a constant (i.e. depends on neither time nor position), 
the continuity equation reduces to: 


V • v = 0 p=const. 

Recall that V • v represents the rate of expansion of fluid elements. “V • v = 0” means that any flow into 
a fluid element is matched by an equal flow out of the fluid element: accumulation of fluid in side any 
volume is negligible small. 


Reynolds Transport Theorem 

In the derivation above, the boundaries of the system were fixed spatial points. Sometimes it is 
convenient to choose a system whose boundaries move. Then the accumulation term in the balance will 


* If the domain V were not arbitrary, we would not be able to say the integrand vanishes for every point 
in the domain. For example: 



sinG dO = 0 



- sin0)r/0 = 0 


does not imply that cos0 = sin0 since the integral vanishes over certain domains, but not all domains. 
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involve time derivatives of volume integrals whose lim its change with time. Similar to Liebnitz mle for 
differentiating an integral whose lim its depend on the differentiation variable, it turns out that:* 


f 




d_ 

dt 


\S{r,t)dV 

K V(t) 


J ^-dV + Js(r,f)(n-w)rfa 

V{t) dt A{t) 


(6) 


where w is the local velocity of the boundary and S(t ) is a tensor of any rank. If w = 0 at all points on 
the boundary, the boundary is stationary and this equation reduces to that employed in our derivation of 
the continuity equation. In the special case in which w equals the local fluid velocity v, this relation is 
called the Reynolds Transport Theorem .♦ 

EXAMPLE: rederive the continuity equation using a control volume whose 
boundaries move with the velocity of the fluid. 

Solution-. If the boundaries of the system move with the same velocity as 
local fluid elements, then fluid elements near the boundary can never cross it 
since the boundary moves with them. Since fluid is not crossing the 
boundary, the system is closed. For a closed system, conservation of mass 
requires: 



or 


_ri_Jmass of] _ q 
dt 1 system J “ 


clM 

dt 


J pdV = 0 

V(t 


(7) 


Notice that we now have to differentiate a volume integral whose lim its of integration depend on the 
variable with respect to which we are differentiating. Applying (6) with w=v (i.e. applying the Reynolds 
Transport Theorem): 


* For a proof, see G: 163-4. 

* Osborne Reynolds (1842-1912), Engineer, bom in Belfast, Northern Ireland, UK. Best known for 
his work in hydrodynamics and hydraulics, he greatly improved centrifugal pumps. The Reynolds 
number takes its name from him. 

* When we say “closed,” we mean no net mass enters or leaves the system; individual molecules might 
cross the boundary as a result of Brownian motion. However, in the absence of concentration 
gradients, as many molecules enter the system by Brownian motion as leave it by Brownian motion, v is 
the mass-averaged velocity. 
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( 




d 

dt 


[ p(r, t)dV 
V(0 


J 


\ ^-dV + f p(r,?)(n-v)<:/<:/ 
J at J 

V(t) A(t) 


which must vanish by (7). Applying the divergence theorem, we can convert the surface integral into a 
volume integral. Combining the two volume integrals, we have 


1 

v{t) 


dp 

dt 


+ V-(pv) 


dV = 0 


which is the same as we had in the previous derivation, except that V is a function of time. However, 
making this hold for all time and all initial V is really the same as holding for all V. The rest of the 
derivation is the same as before. 


Stokes Theorem 

Let v be any (continuously differentiable) vector 
field and choose A to be any (piecewise smooth, 
orientable) open surface. Then 

Jn • ( V x v) da = j> v -dr 

A C 



where C is the closed curve forming the edge of A 

(has direction) and n is the unit normal to A whose sense is related to the direction of C by the “right- 
hand rule”. The above equation is called Stokes Theorem * 


Velocity Circulation: Physical Meaning 

The contour integral appearing in Stokes’ Theorem is an important quantity called velocity 
circulation . We will encounter this quantity in a few lectures when we discuss Kelvin’s Theorem. For 
now, I’d like to use Stokes Theorem to provide some physical meaning to velocity circulation. Using 
Stokes Theorem and the Mean Value Theorem, we can write the following: 


* Sir George Gabriel Stokes (1819-1903): British (Irish bom) mathematician and physicist, known for 
his study of hydrodynamics. Lucasian professor of mathematics at Cambridge University 1849-1903 
(longest-serving Lucasian professor); president of Royal Society (1885-1890). 
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Stokes' Mean Value 

Theorem Theorem 

JLfc/r = |n-(Vx \)da = ^n-(Vx v)^ A = 2(p. n ) A 
C A 

Finally, we note from the meaning of curl that Vxv is twice the angular velocity of fluid elements, so that 
n- Vxv is the normal component of the angular velocity (i.e. normal to the surface A). Thus velocity 
circulation is twice the average angular speed of fluid elements times the area of the surface whose edge 
is the closed contour C. 

Example: Compare “velocity circulation” and “angular 
momentum” for a thin circular disk of fluid undergoing 
solid-body rotation about its axis. 

Solution: Choosing cylindrical coordinates with the z- 
axis aligned with axis of rotation. Solid-body rotation 
corresponds to the following velocity profile (see page 
22 ): 

v = rfleg 

and Vxv = 2f2e- 

Finally the unit normal to the disk surface is n = e_. Then the velocity-circulation integral becomes 

| v -dr = Jn-(Vx v)<r/n = Je- -(2 tie.) da = 2 QjiR 2 

C A A 

According to L&L Vol I* page 25, the angular momentum L of a mass m undergoing motion at velocity 
v is the lever arm r times the linear momentum (p = mv): i.e. L = rxp. Summing this over differential 
fluid mass in our disk with dm = p dV, the net angular momentum of the disk is: 

L = J (r x v)pr/V = pAzj* (r x \)da 
V A 

Since the disk is of uniform thickness Az and density p, we can write the second equation above. If the 
disk is sufficiently thin that we can neglect the z contribution to the position vector, then we can 
approximate r = re r in cylindrical coordinates. 4 Substituting into the second integral above 



* Landau & Lifshitz, Mechanics and Electrodynamics (Course of Theoretical Physics: Vol. 1), 
Pergamon, 1959. 
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R 

= pAzJ [r~D.e z )da = pAzf2e, J r~ 2nrdr = — i? 4 pAz£2e- 


Dividing this by the velocity circulation integral: 

71 R 4 pAz£l 


L z _ 2 


<>vrfr 2 QjiR 2 ^ 


4^=^^p=£- 


V 


4n 


c 


where M is the mass of fluid in the disk. This could be rewritten as 


C 


<\-dr = 4 k — — 
M 


So the velocity-circulation integral is just proportional to the angular momentum per unit mass. 


Derivable from a Scalar Potential 

A very special class of vector fields consists of those vectors for which a scalar field exists such that 
the vector can be represented as the gradient of the scalar: 

Suppose: v = v(r) and/=/(r) 

If/ exists such that: v = V/' 

for all r in some domain, then/(r) is called the scalar potential of v and v is said to be derivable from 
a potential in that domain. 

An example of a vector field which is “derivable from a potential” is the 
gravitational force near sea level: 

Fgrav — - 37yk 

and the associated potential energy is: 

<Kz) = Mgz 



♦ Actually this assumption isn’t necessary since any z-component of r will produce an r-component in 
the cross-product and this r-component will integrate to zero as long as V is a body of rotation about 
the same axis. 
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Note that V<f» = Mgk. 

is identical to the force, except for the sign (introduced by convention). This example also suggests why 
(|) is called the “potential” of v. Not every vector field has a potential. Which do? To answer this, let's 
look at some special properties of such vector fields. 

Property I: if v= V(f) then Vxv=0 (irrotational) 

Proof : Recall that Vx(V(])) = 0 (see HWK #2, Prob. 4e). A vector which has this property is said to be 
irrotational . This name is an allusion to Vxv representing the rotation rate if v is the fluid velocity. 
Vxv=0 means the fluid elements are not rotating. 

Property IP. if v=V(|) then | v -dr = 0 (conservative) 

c 


for any closed contour in the region. 

Proof : Using Property I, we know that Vxv=0. Then we can deduce the value of this closed-contour 
integral from Stokes’ Theorem: 


J-t- dr = | n- da = 0 

C A o 

A vector field which has this property is said to be conservative. This name is an allusion to the special 
case in which v represents a force, like gravity. Then v • dr (force times displacement) represents the 
work required to move the object through the force field. Saying that the contour integral vanishes 
means that the work required to lift a weight can be recovered when the weight falls. In other words, 
energy is conserved. 

If C is open, v=V(|) is still quite useful: 

Property III: let C Q be an open contour connecting points A and B 
Ifv=V(|) then \ Co \ • dr = ()(r s )-(|)(r 4 ) 
for any contour connecting A and B. 

Proof : Note that V(|) • dr = d§ (from our definition of gradient). Th 

| \ dr= | </<|> = <|>(r B )-<|>(r A 

C 0 C 0 

We call this property path independence . Of course, Property II is just a special case of this for 
which A =B so that () ( rg) - <|)(r A ) = 0. 
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Theorem III 

We have just shown that properties I and II are implied by v = Vcf); it turns out that the converse is 
also true, although I’m not going to prove it here. We can distill these properties and their converse into 
a single statement: 


Vx v= 0 
for all r 

► <=> < 

(f> = 4>(r) exists 
such that v = Vcf) 

> <=> < 

| v -dr = 0 for every 
C 

in Region 


in Region 


closed C in Region 


Transpose of a Tensor, Identity Tensor 

The transpose of a tensor x is denoted ^ and is defined so that: 

v • X = • v 

and x • v = v • x f 

for all vectors v. For example: 
if x = ab 

then = ba 

More generally, in terms of scalar components of x, we can write the relationship between a tensor and 
its transpose as: 


Symmetric Tensor. x f = x 

An example of a symmetric tensor is the dyad aa. 

Identity Tensor. Also known as the Idem Factor. Denoted as I and defined so that: 

v ■ I = v = I- v 


for any vector v. Clearly I is symmetric, but in addition, dotting it with another vector gives that vector 
back (like multiplying by one). In any coordinate system, I can be calculated from: 


1 = Vr = — 
= dr 
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where r is the position vector expressed in terms of the unit vectors in that coordinate system. Recalling 
from p6 that gradient can be thought of as the partial derivative with respect to position, I can be 
thought of as the derivative of the position vector with respect to itself. In R.C.C.S., recall that: 

r = xi + yj + zk 

and Vr = LL ^ e / e / 

i j ' 1 

where r,- is the j th component of the position vector r and is the / lh coordinate. In Cartesian 
coordinates, the position vector components are related to the coordinates according to: 

r \~ x \ = r 2 = x 2 = y-> r 3 = x 3 = Z'- 


then 


dr ; 


U 


which is 0 if vtj or 1 if i=j. This leaves: 

Vr = IE*<V 

i j 

so I = ii + |j + kk 

As a partial proof that I has the desired properties which make it the identity tensor, consider dotting it 
with an arbitrary vector v: 

v • (ii+jj+kk) = vii + v- jj + v-kk 
= (v • i)i + (v • j)j + (v • k)k 
= v x i + Vyj + v z k = v 

Thus we have shown that v • I=v, as advertised. 


Divergence of a Tensor 

In presenting the corollaries to the Divergence Theorem, we have already introduced the divergence 
of a tensor. This quantity is defined just like divergence of a vector. 


* This expression for the identity tensor is valid for any set of orthonormal unit vectors (not just the 
Cartesian ones for which we have derived it here). 
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V-x 


lim 
V ->0 


V 


Jn -Tela > 
A 


Note that this definition uses a pre-dot not 
post-dot. In R.C.C.S. 

I = SEXy-e/ey 

On the x=x 0 face: 

n = -e i 


n-x 


I 


E T y e l' e i e j 

j 


“E 



Similarly, on the x=x 0 +Ax face, we obtain: 

n = -t-e i 


"i = E 

j 


x a +At 


After integrating over the area, we obtain: 


J n-x da = £ {x i j \ X/j +Ax -x | j \ Xo }e , AyAz 
Aj + A 2 7 


Dividing by V: 




Ai + A 2 j 


Ar 


Taking the lim it as V^0: 


lim 

V->0 


1 r 

— n- xda 
V . J . = 


A\ + A 2 


3xi j 
' Oa Cj 


Adding on similar contributions from the y=const and z=const faces: 
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V-T 


IE 

i j 


1 


V-T 


^ z xx | ^ X yx | 
dx dy 


^tx 

dz 


\ 

i + 


dz dz 


xy 


+ 


dz 3 

XV + ° L zy 


dx dy dz 


j + 


J 


dVvz _j_ d^yz _j_ d^zz 
dx dy dz 


\ 

k 


Introduction to Continuum Mechanics* 

In this course, we will treat fluids like continua; in other words, we are going to ignore the molecular 
granularity of matter. This is an assumption which we, as engineers, often make in describing transport 
of heat, mass, or momentum although we don’t always state this assumption explicitly. To make the 
nature of this assumption clearer, it might help to discuss the alternative. 

Fluids are composed of molecules. In principle, if you tell me the initial location of every molecule in 
the fluid and its initial velocity, I can compute the position and velocity at some later time using 
Newton’s laws of motion (i.e. F = ma). The difficulty with this approach is that the number of 
molecules in any volume of fluid of interest to us make such a detailed calculation impractical. For 
example: 


1 cm3 0 f water — > 3.3xl0 22 molecules — > 10 million years 

Even with a computer operating at 100 mfops, it would take 10 mi lli on years to do just one 
multiplication for each molecule. Molecules of a liquid collide on the average of once every 10 -12 
seconds. To describe one second of real behavior, I would need 10 12 x 10 mi lli on years. Clearly, this 
is an absurd length of time. Although computers get faster every year, this will remain an absurdly long 
time for the foreseeable future. The alternative is: 


Continuum Hypothesis 

A detailed description at the molecular level is not required in order to predict macroscopic 
behavior of any material. For example, it is not necessary to know the precise location of every 
molecule of fluid; it turns out that all that is needed for most applications is the distribution of mass 
described by the density profile p(r) of molecules in some region: 


* Reference: G.E. Mase, "Continuum Mechanics," Schaum's Outline Series, McGraw-Hill, 1970. 
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p(-> = ;: 0 i>] 

where m ( - is the mass of molecule i, the sum is over all the molecules 
inside surface A, and V is the volume. In fluid mechanics, as in heat and 
mass transfer, we make an assumption known as the "continuum 



hypothesis." 

Basically this assumption is that the lim it above will 
converge long before the dimensions of V shrink to 
molecular size. 

Similarly, we don't need to know the translational, 
rotational, vibrational and electronic energy of each 
molecule. We usually need only to know the internal 
energy per unit volume as a function of position, which in 
turn, manifests itself macroscopically as temperature. 




Q 1 1 1 — 

|0rll> 


-f ► V (cm 3 ) 

10 ° 


A more precise statement of the continuum hypothesis is: 


Continuum Hypothesis - the region to be described can be subdivided into a set of 
(infimtesimal) volume elements, each of which simultaneously: 


1) is small enough to be considered uniform (i.e. any spatial variations in properties — such as p, v, 
T, p — inside the volume element are negligible); and 


2) is large enough to contain a statistically large number of molecules. 


In other words, we are assuming that dV exists such that the two conditions above are both satisfied. 
Materials which obey this “hypothesis” as said to behave as a continuum. Generally, the continuum 
hypothesis works well provided all the dimensions of the system are large compared to molecular size. 
An example of a situation in which the continuum hypothesis does not work is the flow of dilute gases in 
small pores, where the mean free path (for the collision of molecules) is comparable to the dimensions of 
the pore. This situation is called “Knudsen diffusion.” 


The basic problem in continuum mechanics is to describe the response of material to stress. A 
quantitative statement of that response is known as: 


Constitutive Equation - model which describes how a material will respond to stress. 
Familiar examples of constitutive equations include: 
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1) Hooke's law of elasticity (solids) 

2) Newton's law of viscosity (fluids) 

Many materials, like toothpaste and polymer melts, have characteristics of both solids and fluids and do 
not obey either of these simple "laws." Such fluids are called viscoelastic. 


Classification of Forces 

Having derived an equation (the Continuity equation) to describe the relationship among the 
variables which is imposed by conservation of mass, the remaining fundamental principle of physics is 
Newton’s second law (E ( -F, = ma) which, as it turns out, is equivalent to conservation of momentum. 
To apply this principle, we will need to list the forces which can act on fluid systems. Forces tend to fall 
into one of two different categories, depending on the range over which they act: long-range (compared 
to molecular size) forces are computed as volume integrals (called “body forces”) and short-range 
forces are computed as surface integrals (called “surface forces”). 

Of course, gravitational forces have the longest range of any known force. For example, gravitational 
forces between planets and the sun determine their orbits. In particular, all fluid elements (not just those 
at the system boundary) feel a gravitational force of interaction with the rest of the universe outside the 
system boundaries. Thus gravity is a “body force.” 

body forces : those which act on every fluid element in body (e.g. gravity): 

dF g = (dm ) g = pg dV 

At the other end of the spectrum are forces which 
have very short range. If the range is of molecular 
dimensions, then only fluid elements experience a 
nonzero interaction with the universe outside the 
system. Although interior fluid elements might 
interact with one another through this short-range 
force, this interaction is not considered in a force 
balance, because the “action” and “reaction” 
forces cancel, leaving no net contribution to the 
force on the system. When only surface elements 
feel a particular force from outside, that force is called a “surface force.” 

At the molecular scale, pressure arises from the momentum transferred during collisions between 
molecules outside and molecules in side the system. Since only surface molecules will be struck from 
outside, pressure is a surface force. 


short-range interactions heLween 
’iimefi.Oi. 1 ’' elements cancel 



short-range interactions between 
"surface" ckmcnLs do nol cancel 
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surface forces : those which act only on surfaces 
(including mathematical boundaries) 

One example is hydrostatic pressure : 

dF p = -pn da 

where dF p is the force exerted on the system (through 
the surface element da) by the fluid outside, and n is a 
unit outward (to system) normal. For a proof that this 
is the correct form for hydrostatic pressure, see 
Batchelor, Section 1.3. 


n 



Hydrostatic Equilibrium 

If our material is a fluid and if it is at rest (no velocity and no acceleration), then gravity and 
hydrostatic pressure forces are usually the only forces acting on the system. At equilibrium, the forces 
must be balanced. Thus Newton’s* 2 nd law, which generally requires ^ F, = A/a . reduces to 

i 

y F, =0 at mechanical equilibrium. 

i 


In our case, this means 


F + F =0 

g p 


F g = \vP%dV 
F„ = -\ A npda = - \yVpdV 

To obtain this last result, we applied one of the corollaries of the Divergence Theorem. Substituting 
back into the force balance and combining the two volume integrals leads to: 


F + F 

r g r p 


Iy[pg-V//|dV = 0 
Since V is arbitrary, we conclude that the integrand vanishes: 

Vp = pg 


air 


water 


* Sir Isaac Newton (1642-1727), English mathematician and nai 

considered by many the greatest scientist of all time; invented differentia 

theories of universal gravitation, terrestrial mechanics, and color. 


lower p 

Vp 




higher p 


it); 

.he 
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This says that the pressure increases in the direction of the acceleration of gravity (downward), which 
correctly describes (for example) how the pressure increases with depth in an ocean. 


Flow of Ideal Fluids 


Now let's consider fluids in motion. The simplest analysis is for: 

ideal fluid - deformation of fluid elements is an isentropic process (i.e. adiabatic and 
reversible): 


(l = 0 and k = 0 

where p is the viscosity and k is the thermal conductivity. Generally this means that any viscous forces 
are negligible (since viscous forced represent friction arising between fluid elements and friction gives 
rise to irreversibility). Furthermore, to keep the process adiabatic, the thermal conductivity must also be 
negligible. 


Euler's Equation 

Suppose these conditions on the fluid are met. Thus consider the isentropic deformation of an ideal 
fluid for an arbitrary macroscopic system. In addition to pressure and gravity, we must also consider 
inertia when the system accelerates. Newton's law requires: 

Ma = E ; F ; 

Let if/) denote the trajectory of one particular fluid 
element inside the system. Then the velocity of the 
fluid element is: 

Dr 

v = 

Dt 

while the acceleration is: 

Dv 

a = 

Dt 


trajectory 



We use the material derivative here, since we are following a particular material point. Multiplying the 
acceleration by the mass of the fluid element gives the inertia: 


(r/m)a = p (dV) 


Dv 

~Dt 
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To get the net inertia of the entire system, we must repeat this calculation for each of the fluid elements 
composing the system and add them up: 


(Ma) = | p 


V 


Dv 

Dt 


dV 


I’ve put carets (<...>) around the Ma to indicate that this is the position-average inertia of the system 
(since the local fc/m)a varies from point to point within a general system). Newton's second law 
requires us to equate this with the net force acting on the system: 


JP 


Dv 

Dt 


dV = F 


p 


+ F 


Jpgr/V-J n pda 
12 3 


v 


jVpdV 

V 


Using the divergence theorem to convert the surface integral into a volume integral, we have three 
volume integrals over the same domain. Combining these three volume integrals leaves: 


r Dv ^ 

J P — -Pg + Vp 


\dV = 0 


Since this must hold for any choice of V, the integrand must vanish at each point in the domain. After 
dividing by p : 


Dv 1 

— =g Vp 

Dt p 


(8) 


which is called Eider's Equation (1755).* 

Significance: When combined with a statement of continuity, Euler’s equation of 
motion provides as many equations as unknowns. 


Another relationship among the unknowns is the continuity equation (see page 27), which comes 
from the mass balance. 


* Euler, Leonhard, 1707-83, Swiss mathematician. The most prolific mathematician who ever lived, he 
worked at the St. Petersburg Academy of Sciences, Russia (1727-41, 1766-83), and at the Berlin 
Academy (1741-66). He contributed to areas of both pure and applied mathematics, including 
calculus, analysis, number theory, topology, algebra, geometry, trigonometry, analytical mechanics, 
hydrodynamics, and the theory of the moon's motion. 
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¥ + V'(p4 = o 

For an incompressible fluid, p = const. w.r.t. both time and position. Then the continuity equation 
reduces to: 

imcompressible fluid: V • v = 0, (9) 

To see that we now have as many equations as unknowns, note that the unknowns in (8) and (9) 
are 

unknowns: 4 scalars v and p 

which represents the 3 scalar components of v plus p, for a total of 4 scalar unknowns. To evaluate 
these unknowns, we have equations (8) and (9): 

equations: 4 scalars Euler + continuity 

but Euler’s equation (8) is a vector equation, which can be expanded into 3 independent scalar 
equations. When added to continuity (a scalar equation), we obtain a total of 4 independent scalar 
equations, the same number as of scalar unknowns. Thus we are now in position to begin solving 
problems involving fluid flow. We will call (8) and (9) “Euler’s equations of motion for incompressible 
fluids.” 


EXAMPLE. Water in a partially fill ed tank undergoes 
uniform* acceleration a in the horizontal plane. Find the 
angle 0 of inclination of the water’s surface with respect to 
the horizontal plane. 



a 

► 


2^2 


Solution. The key to solving this problem is to recognize that, 
regardless of the angle of inclination, the pressure is equal to 
1 atm everywhere on the free surface. Then the pressure 
gradient Vp must be normal to this plane. If we can find the 

orientation of Vp we will have the orientation of the free surface. Recall Euler's equation of motion for 
an ideal fluid: 


D\ 1 w 

— = g --Vp 
Dt p 


* By “uniform acceleration” I mean the same acceleration is experienced at each point and at each time. 
For a fluid, uniformity at each position occurs only in the steady state after a transient which is 
nonuniform. 
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D\/Dt is just the acceleration of the fluid in a stationary reference frame. At 
steady state, all of the fluid will undergo the same uniform acceleration as the 
tank; so D\!Dt is just a. Solving for the gradient, we have 

-Vp = g- a 

P 


const 

r surface 


Xp 



30 


Using vector addition in the drawing at right, we can see that the angle of 
inclination of the free surface (relative to the horizon) is just 


0 = tan 


-1 


r a ' 


We were lucky in the previous example, because we knew the left-hand side of (8), so instead of 4 
scalar unknowns, we only had one: p. The solution was relatively easy. In the more general problem, 
the left-hand side of (8) is an unknown nonlinear partial differential equation: 

dv 1 

— + v-Vv = g Vp (10) 

dt p 

In this form, we have expanded Dv/Dt using the relationship between material derivative and partial 
derivative (see page 15). Now we have 4 scalar unknowns: the three scalar components of v and 
pressure: v x , v y , v, and p. Coupled with the continuity equation (for an incompressible fluid) 

V- v = 0 

Euler’s equation also gives us 4 scalar equations. One important class of solutions has the form v = Vp, 
which is called “potential flow.” In the next section, we discuss how this form comes about and identify 
which physical problems have this form. 


Kelvin's Theorem 

An important precursor to the theory of potential flow is the principle of conservation of 
circulation . Before stating this principle, let me define a quantity which Landau & Lifshitz* call the 
velocity circulation : 


* L.D. Landau and E.M. Lifshitz, Fluid Mechanics (Vol. 6 of a Course of Theoretical Physics), 
Pergamon, New York, 1959. 
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T = j> vr/r 

C 

for any closed contour. Recall that we showed on page 
30 that this contour integral is associated with the average 
angular momentum of fluid elements located on the surface 
whose edge is C. Kelvin 4 showed that this velocity 
circulation is conserved: 



contours 

of material points 
at different times 



C(t 2 ) 


trajectory of a 
material point 

a material point 


for any set of material points forming a closed contour in 
an ideal fluid. This result is called Kelvin 's Theorem. 


Partial proof : y Since the contour C is composed of material points, the time derivative of this contour 
integral is like the material derivative: 


dT_ 

dt 


Dr • r Dv 
— ±jL-<7r = -L = — dr 
Dt ~ ~n t 


C 


C 


Dt 


Since we are always integrating over the same set of material points, a boundary term does not arise 
when we interchange integration and differentiation operators, although the set of spatial points is time- 
dependent: C = C(t). Of course, we have not rigorously shown this step to be valid, thus we only claim 
the proof is “partial.” Next, we substitute Euler’s equation and write each term as the gradient: 


dT 

dt 



i 3 f 

— vp -dr = i 

P J J 


'-v^-v4* = -fv|^ + k| 


C 


P) 


c 


dr = 0 


In the second equality above, we have introduced the potential energy per unit mass, (|) g . Recall that 
gravity is a conservative force field (see page 32). For an object of constant mass (e.g- a brick), 
Theorem IH guarantees that a scalar field <f> (r) exists such that F g = mg = -V<]> . For an object of 


♦ Ford Kelvin (William Thomson), 1st Baron, 1824-1907, British mathematician and physicist; b. 
Ireland. He was professor (1846-99) of natural philosophy at the Univ. of Glasgow. His work in 
thermodynamics coordinating the various existing theories of heat established the law of the conservation 
of energy as proposed by James Joule. He discovered what is now called the “Thomson effect” in 
thermoelectricity and introduced the Kelvin scale, or absolute scale, of temperature. His work on the 
transmission of messages by undersea cables made him a leading authority in this field. 

v For a more rigorous proof, see Batchelor p269. For a more intuitive proof, see F&F, pl5. 
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constant mass, we could divide both sides of this equation by m and bring m inside the gradient 
operator: g = -V(<( )/m). Similarly, a differential volume element dV would have a differential mass dm 
and a differential potential energy <:/<]). such that g = - V(d<\>/din). So generally we can write: 

g = -V(f>g (11) 

, , potential energy 

where cp „ = . 

5 mass 


The next-to-last expression (in the equation for dT/dt ) must vanish, because it has the form V.s’ • dr = ds : 
integrating any total differential around a closed contour yields zero. Thus 

v -dr = const, w.r.t. time 

C 

Keep in mind that this applies only if C is composed of material points and only for ideal flow. Since C 
is composed of material points (which in general move with different velocities), the contour may change 
shape or move. Given the meaning of this contour integral (T = angular momentum), this result implies 
that (in the absence of friction) angular momentum of fluid is a constant. In general, we change the 
angular momentum of some object by applying a torque. So this result (i.e. Kelvin’s theorem) means 
that (in the absence of friction) there is no way to apply a torque to ideal fluid 
elements. 

If you think about it, this makes sense: the usual way to apply a torque (with our 
hands to a cylinder, say) is to hold the cylinder between our hands and then move 
our hands in opposite directions, as shown in the sketch at right. We thus rely on 
friction between our hands and the cylinder to exert the torque. If the cylinder 
were greased and our hands slipped over its surface, we would not be able to 
apply the torque. This is the essence of what Kelvin’s theorem is saying. 



IRROTATIONAL FLOW OF AN INCOMPRESSIBLE FLUID 


As an example, consider towing a submerged object through an “ideal fluid” which is otherwise 
stagnant. 


v(r,f=0) = 0 



submerged object 

mov ing Thro ugh 

stagnant fluid 


Consider an arbitrary closed contour in the fluid far from 

the disturbance caused by the motion of the submerged object. The contour integral vanishes since v 
vanishes: 


| v -<r/r = 0 for 1=0 since v = 0 

C 
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for every such contour C. Now at some later time the submerged object moves into the vicinity of C 
which causes v to be nonzero. Despite this, Kelvin’s Theorem still requires 


| v -dr = 0 for all l although v # 0 

C 


This is also hue for every closed contour in the region (since every contour initially had zero value for 
this integral). 

Applying Theorem III: Vxv = 0 for all rj (12) 

which is sometimes called the persistence of irrotationality . Also from Theorem HI, we know (fir,/ 1 ) 
exists such that: 


v = V<|) 


(13) 


where (f> is called the velocity potential. 


Significance: Knowing that the solution has the form given by (13) allows us to 
decouple the four scalar equations represented by (8) and (9): 


Substituting (13) into (9) yields Laplace’s equation in the velocity potential: 

(13) into (9): V 2 p = 0 

Instead of 4 equations in 4 unknowns, we now have a single equation which can be solved for cf> and v 
= V(|), without any coupling to Euler’s equation (8). Although Euler was the first to suggest this 
approach, this is called Laplace's Equation* after the French mathematician who solved this equation 
in so many cases. 

Knowing the velocity profile v, we can now determine the pressure profile p from Euler's equation: 

^L + vVv = g--Vp (14) 

dt p 

We will now integrate this vector equation to obtain a single scalar equation for the pressure profile. 
Each term in (14) can be expressed as a gradient of something. For example, we’ve already seen in 
(11) that: 


* Pierre Simon de Laplace (1749-1827), French mathematician and astronomer, noted for his theory of 
a nebular origin of the solar system and his investigations into gravity and the stability of planetary 
motion. 
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♦ “Identity A.3” in this equation refers to one of the mathematical identities summarized on the handout 
titled “Useful Identities in Vector Notation”. 

♦Daniel Bernoulli (1700-1782, Swiss) has often been called the first mathematical physicist; the teacher 
of Leonhard Euler. His greatest work was his Hydrodynamica (1738), which included the principle 
now known as Bernoulli's principle and anticipated the law of conservation of energy and the kinetic- 
molecular theory of gases developed a century later. He also made important contributions to 
probability theory, astronomy, and the theory of differential equations. 
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Potential Flow Around a Sphere 

To solve a typical problem involving potential flow, we would first solve Laplace's equation to 
obtain the velocity profile and then we can evaluate the pressure profile using Bernoulli's equation. Let's 
illustrate the procedure using an example: 

EXAMPLE: Find the velocity and pressure profiles for potential flow 
caused by a sphere of radius R moving through a stagnant fluid with v — ® 
velocity U. 

Solution: If the fluid behaves ideally, it undergoes potential flow and 
the velocity profile must satisfy Laplace’s equation: 

P.D.E. : V 2 4> = 0 

Boundary conditions can be formulated by recognizing that fluid far from the sphere is unperturbed: 

b.c. #1 : v = 0 far from sphere 

while fluid near the sphere cannot penetrate the sphere. To express this mathematically, recall our 
“bucket-and-stopwatch” method for defining fluid velocity (see page 13). Modifying it slightly to 
account for the movement of the surface element at velocity U, the flowrate across a surface element of 
area da is given by: 



dq = n-(v -U )da = 0 


For an impenetrable sphere, the flowrate must vanish 

b.c. #2 : n • (v - U) =0 on sphere 


To solve this problem, we adopt a new reference frame 
in which the origin moves along with the center of the 
sphere. It turns out that the PDE does not change upon 
this shift in reference frame for velocity: the new velocity 
potential must also satisfy Laplace's equation: 

P.D.E. : V 2 4> = 0 

However, the boundary conditions are changed. 

In this moving coordinate system, the sphere — * 

appears to be stationary and the fluid at infinity is » 

undergoing uniform flow 

b.c. #1 : v — > -U = Uk as r— > °° (19) * 


Let reference frame move with object: 




uniform flow 
around 

yiiLlionary object 
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where U is the velocity of the sphere in the original stationary reference frame. 

For a stationary sphere, no fluid entering the sphere means dq= n • v da = 0, or 
b.c. #2 : n- v = v r = 0 at r=R 

Next we rewrite the b.c.’s in terms of the velocity potential in spherical coordinates. In terms of velocity 
potential, (19) becomes: 

V9 = Uk 

Now we need to translate k into the unit vectors in spherical coordinates.* 



Referring to the figure above (see page 8), we note that e r e e , and e. = k all he in the same <t>=const 
plane (shaded region of left-hand figure above). If we shift all three unit vectors to the origin (recall that 
the origin is not part of the definition of any vector) and re-orient the <I>=const plane to coincide with the 
plane of the page, then we get the figure above at right, from trigonometry of the right triangle it is 
apparent that 

k =(cos0)e,. -(sin0)ee 

Thus the b.c. can be written as 

as r— > oo; V<)) — > ^ 2o§0e ; . l4 / 2'i60 e e 

3<j) 1 3(j) 

dr r 30 

Equating corresponding components: 

^- = U cos 0 and = U sin0 

or r o0 


* Here <E> is the spherical coordinate, while (|) is the velocity potential. For more on spherical 
coordinates, see BS&L p733. 
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Integrating either PDE leads to 

b.c. #1 : 4> — » I/rcosO + const as <» 

where we have arbitrarily selected a value of zero for this “const”.* In order to translate b.c. #2, use 
tables (see spherical coords 4 ) to express the gradient in spherical coordinates: 


v =V(|) 

a<f> i a* i a<b 

v r e r + v e % + v <j> e <t> - e r + e e + — . e ft » 

or rod rsinOckP 


Dotting both sides by the unit vector n = e r , then using b.c. #2: 


b.c. #2: 


v r = 3 4>/Dr = 0 at r=R 


We look for a solution which is independent of <f>: 

<|> = <Kc0) 

This imphes that the fluid does not have any d>- 
component of velocity. In other words, the trajectory 
of any fluid element remains entirely on a single 
<f>=const surface. The sketch at right shows some 
edge views of <J>=const planes (looking along the z-axis 
which passes through the center of the sphere). 

From p740 of BS&L, we have V 2 (|) in spherical 
coordinates: 



edge view of 
<I> - const 
surfaces 


end view of 
line parallel 
i -/-axis 


j__a/ 

r 2 3r v 


a<l> 

dr 


+ - 


1 


sin0 


r _ sin0 ^0 v 
d§ldr = 0 at r=R 
<j> = Urc os0 as r— > oo 


3(f) 

30 


= 0 


( 20 ) 

( 21 ) 

( 22 ) 


* Lik e most energies in thermodynamics, the reference state for potential is arbitrary and can be chosen 
solely for convenience. 

♦ http://www.andrew.cmu.edu/course/06-703W ops_sph.pdf 
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Although there are systematic procedures, like “separation of variables,” for solving P.D.E.'s which 
work in many cases, one should first check to see if there is a simple solution. The problem appears 
formidible, but notice that the trivial solution: 


<|> = o? 

satisfies the P.D.E. and the b.c. at r=R. Unfortunately, it does not satisfy the b.c. at r — > °°, so we will 
have to try another guess. Since the failure occurred with the b.c. at r—> oo, we look there for the form 
of our second guess: 


4> = Arc os0 ? 

This too satisfies the P.D.E. and the b.c. at r— > °° (provided A=U) but to satisfy the b.c. at r=R, A must 
be chosen as 0. The tells us that A should have different values at different r's. So we try a third guess 
which is slightly more general than the second: 

<Kr,0) =/(r)cos0 (23) 

Substituting (23) into (20)- (22), we find that cosO cancels out, leaving: 

r 2 f" + 2rf'-2f=0 
f ' = 0 at r=R 
f= Ur as r—> oo 

Thus we have reduced the problem of solving the P.D.E. to one of solving an O.D.E. We recognize this 
O.D.E. to be a Cauchy-Euler equation,* which always has at least one solution of the form /= r n . The 
general solution turns out to be: 


f(r) = A r- + Br 


* The general form of an -order Cauchy-Euler equation is 


d y 

L a„x n — — = 0 . At least one of the 

rlr n 

n = 0 


N linearly independent solutions has the form y(x) = Ax a . Substituting this form leads to 


N 

a 0 Ax a + £\/ ;) a(a - l)(a - 2)K (a - n + l)A.r a = 0 

n — 1 


Dividing out the common factor Ax a , we obtain an M h -order polynomial for a. Each distinct root of 
the polynomial leads to a separate solution. In this example, the roots are a = -2 and +1 . 
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Using the b.c.'s, we can evaluate A and B. The particular solution to this problem is: 
f(r) = U 

Substituting this back into (23) leaves: 


' t^ 
r+ 2- 
V r ) 




x R 3 ^ 
r+ 2- 
V r J 


3d) 

v r =/ = U 
or 


1- 


' R_' 

\r) 


COS0 


COS0 


1 ^ TT 
V ° = ^ = ~ U 


f i R 3 ^ 


1 + 


2 3 


sinO 


if 


v r = 0 


+ 


v 0 (r,e = f) 


i r 

f r*- 



v r (r,0) 

V0=O 


v e (r,9 = f) 


Notice that for 0 =tc (0=0), v e =0 but v r decreases (in absolute 

magnitude) from -U at r=°° to 0 at r=R. At 0 =tc/ 2, v r =0 but v e increases from U to (3/2)67. This 
increase is necessary to make up for the decrease in flow caused by the sphere blocking part of the flow 
path. 

Having solved for the velocity profile, we can determine the 
pressure profile from Bernoulli's equation (18). Assuming steady 
state: 

2 

V p 

1 1- <p „ = const 

2 p 



independent of position. The "const" can be evaluated using any point where we know both the velocity 
and the pressure. Suppose the pressure of the undisturbed fluid is known in the reference plane for the 
gravitational potential (() „ = 0): 


at r— > °o, p = 0: 

p = Poe and v 2 = U 2 

Thus 

V 2 P , U 2 Poo 

— + — + <b. = + 

2 p 2 p 

or 

P (r,e) = Poo -p4> ? +^p (u 2 -v 2 ) 

where 

V 2 =yv=[v,.(r,e)]- +[v 0 (r,e)]~ 
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Substituting the known velocity profile, we obtain the pressure profile. Let's focus on the pressure 
profile over the surface of the sphere: 

forr=R. p(R,0) = +Ip 

hydrostatic head dynamic head 

Ph{ R ^) Pd{ R ,Q) 

In the sketch at right, we plot the dynamic pressure (dropping the 
contribution from hydrostatic equilibrium). Note the location of 
regions having high and low pressure. The sphere is being pushed in 
at the poles and allowed to expand at the equator. This is why a 
large bubble rising through stagnant water tends to become distorted 
from spherical shape. Such bubbles tend to become extended in the 
horizontal plane and compressed in the vertical direction by the 
higher pressure. 



after dcibrnuiLjun 



d'Alembert's Paradox 

What is the net force on a rigid sphere owing the pressure profile developed by potential flow 
around it? The answer turns out to be: 

F p = -jnpda = -|n [p h + p c i)da = -j>n p h da - |n p d da = -pgV (24) 

a a 142 43 142 43 

pgV 0 

where V is the volume of the sphere. 

Proof : first consider the contribution from dynamic pressure: 

Pd =iP^ 2 (9cos 2 0-5) 
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The n in (24) is a unit vector normal to the surface, pointing outward. In spherical coordinates, with 
their origin at the center of the sphere, this vector is the unit vector in the r-direction: 


n = e r (0,(|)) 

which direction depends on location on the surface. Although its length is a constant, its direction varies 
with position of the sphere; thus n cannot be treated as a constant. Anticipating that any net force will 
be parallel to the direction of fluid flow, we dot both sides of (24) by k: 

F p = k ■ F p = - J k-n pda 

A cos0 

For the contribution from dynamic pressure, p d (r=R,Q) depends solely 
on 0, so we choose the strip of width R <70 and length 27t/?sin0 as our 
differential area da. On this strip 0 is virtually constant. 

n 

F d P = -Jcos0 (2^7?sin0)(7?r/0) 

0 /(cos0) 

n 

= 2nR 2 Jcos0 /(cos0) (j-^i^0|^) 

0 d cos0 

x=cos0 . 

= 2tt 7?- x±pt/“ )x[9x 2 -5\dx = 0 

1 44 2 4 43 

o 



Hie nonzero contribution from (24) comes from the hydrostatic head. We will leave this calculation as 
an exercise for the reader (HWK #4, Prob. 1). The net force due to pressure is 


F P =~ I (Poo-p<t> g )nda- | ^pf/ 2 (9cos 2 0-5) 


nda =-pgV 


l s Bj3 e 2jF 442 4 4 4 43 f h f e 4 44442444443 

-pgV 0 


The net force on the sphere is the sum of its weight and the net pressure force: 

F p +F s =-pgV + p,gV = (p s -p)gV r *0 


which is the difference in weight of the sphere and the weight of the fluid displaced by it. Because the 
pressure force is independent of the speed of motion of the sphere through the fluid, the particle will 
continue to accelerate forever, without ever reaching a ‘terminal velocity.” Of course, experiments 
show that falling particles reach a terminal velocity which implies that the gravitational force is balanced 
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by some other force. The other force is fluid drag (friction), which is not predicted by potential flow. 
This serious discrepancy between the predictions of potential flow theory and experiment is known as: 

d'Alembert's * paradox - potential flow predicts no drag but experiments indicate drag. 


Despite this, potential flow is still useful: 


Uses of potential flow - predicts lift (but not drag) on 
streamline objects moving through stagnant fluid 
at high Reynolds numbers (but still sub-sonic, i.e. 
v « c). 



• correctly predicts v(r ,t) and p(rj) except very near the surface 

of the object (i.e. inside boundary layer) and in wake [see pressure profile on airfoil shown below]. 


• for asymmetric shapes (e.g. airplane wings), it correcdy predicts a lift. [See lift force shown below as 
a function of the angle of attack.] 



Figs. 1.12 and 1.13 taken from Schlichting, 6th ed., p22f. 


* Jean le Rond d’Alembert (1717-83), French mathematician and philosopher, a leading figure of the 
Enlightenment. His treatise on dynamics (1743) enunciated d'Alembert's principle, which permitted the 
reduction of a problem in dynamics to one in statics. He did important work on the mechanics of rigid 
bodies, the motions of fluids and vibrating strings, and the three-body problem in celestial mechanics. 
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Stream Function 

Before we resolve d’Alembert’s paradox by adding viscous forces, let's step back for a moment 
and review what we have accomplished for potential (or irrotational) flow. The mathematical problem 
might be stated as: find v(x,f) such that: 



Vxv = 0 

(25) 

and 

< 

■< 

il 

o 

(26) 


Eqs. (25) and (26) represent four partial differential equations in the components of the unknown vector 
field * Recognizing that the solution is derivable from a potential allows us compress these four 
equations into one scalar equation in one unknown: 

v=V(f>: v2(J> = 0 


which is quite a remarkable trick. The potential is not the only scalar field which a vector field can be 
expressed in terms of. Velocity can also be expressed in terms of a stream function. 

Potential (<j» ) — a scalar field whose relationship to v is carefully selected to automatically satisfy 
irrotationality 


v = V(f> -> Vxv = 0 

Whereas the relation between velocity and scalar potential is chosen to automatically satisfy (25), the 
relationship between velocity and stream function is chosen to automatically satisfy (26): 


v = f(\|/) — > V • v = 0 


Stream Function (\|/) — a scalar field whose relationship to v is carefully selected to 
automatically satisfy continuity . 

It turns out that it is sufficient to express v as the curl of another vector. According to Identity C.6, the 
divergence of the curl of any vector is zero (See HWK #2, Prob. 4d): 


v = Vxu V- v= 0 


A vector which can be expressed as the curl of another vector is said to be solenoidal. u is called the 
vector potential of v. Of course, knowing that v = Vxu isn't always of much help because we just 


* Although it might appear that we have overspecified the problem by specifying both the divergence 
and curl (which represent four scalar equations in the three components of v), this turns out not to be 
hue. In general, both the divergence and curl must be specified throughout some region in space before 
the vector field can be determined in that region. 
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trade one unknown vector for another. Fortunately, there are several broad classes of flows for which 
the form of the vector potential is known: 


Two-D Flows 


When nothing happens along one of the three directions in R.C.C.S., we have 2-D flow. 

v = v x (x,y)i + v y (x,y)j 

or v z = 0, B/Bz = 0 

identity 

E.3 

v = V x [\\i(x,y)e z ] = + 

For such a flow,* 0 


3\1/ d\|/ 

— +— 1 L e v 

dx x By y 


3\)/ 


xe * = aU^5 

-e„ 


3\i/ 


(27) 


In terms of its scalar components, the velocity is: 


v * = 


3x|/ 

By 


3x)/ 


Bx 


v z =0 


(28) 


Next we substitute this form for v into (26): 

identity 

C.6 

V-v = V-[Vx(\|/e-)] = 0 

which automatically satisfies continuity (26), for any choice of \|/(.v',y). The scalar field x|/(.x,y) is called 
the stream function . For irrotational flow, the problem would be to determine \|/(a',v) such that (25) is 
satisfied: 


Vx[Vx(\(/e-)] = 0 

We can reduce this to a scalar equation. Using identity E.5 from handout: 

Vx[Vx(ye,)] = V[V • (x|/e,)] - V 2 ( ¥ e z ) 


* “Identity E.3” in this equation refers to one of the mathematical identities summarized on the handout 
titled “Useful Identities in Vector Notation”. 
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but V • (\) ie z ) = d\\r/dz = 0 

and v2(\|/e„) = (v2\j/)e_ 

Thus Vxv = -(v2\(/)e_ 

So for irrotational flow, the streamfunction must also satisfy Laplace's equation: 

Vxv = 0 : V 2 \)/ = 0 

U nlik e the scalar potential, the streamfunction can be used in all 2-D flows, including those for which the 
flow is not irrotational. Indeed, we will use the streamfunction to solve Stokes flow of a viscous fluid 
around a sphere, in which the fluid is not even ideal. 


Axisymmetric Flow (Cylindrical) 


Another general class of flows for which a streamfunction exists is axisymmetric flow. In cylindrical 
coordinates (r,0,z), this corresponds to: 


v = v,.(r,z)e r + v z (r,z)e z 

or v e = 0 and 3/90 = 0 

Then V • v = 0 can be satisfied by seeking v of the form: 

v = Vx[/(r,z)e 9 ] 


or 


v = Vx 




The second expression usually leads to somewhat simpler expressions for Vxv and is the one used in 
the table on pl31 of BS&L: 


1 3 \|/ _ 1 3 \)/ 

r dz r dr 


Computing the curl in cylindrical coordinates and setting it equal to zero leads to the following PDE in \(/ 
(the details are left as an exercise): 


V 


x v = 


6 2 3 

E~\\i ] 

V r J 


e fl =0 


E 2 \\f = 0 
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where 


Note that E 2 \\r * : 


_ i d 2 \i/ 1 d\i/ d 2 \i/ 

dr 2 r dr dz 2 


V 2 y(r, z) 


3 2 \(/ 1 3\)/ d 2 \|/ 

dr 2 r dr dz 2 


Axisymmetric Flow (Spherical) 

In spherical coordinates (r,0,())), axisymmetric flow means 

v = v ; .(r,0)e r + v e (r,0)e e 

or = 0 and d/d<\> = 0 

where cf) is the a zim uthal angle. Then V • v=0 can be satisfied by seeking v of the form: 

v = Vx[y '(r,0)e^] 

1 3\i/ 1 3ur 

— — e,. -e 0 

V r v 0 

Again, the second expression is the one used in the table on pl31 of BS&L (except the signs are 
reversed). Taking the curl (HWK #4, Prob. 6a): 


which requires 
where 

Note that E 2 \\t * V^\\i 


V x v = 


rsin0 
E 2 \ r = 0 


* E 2 \\i = 0 


1 d\|/ 

3 r 2 r 2 dotsinOdO, 


—,2 d 2 \|/ sin0 d 

E \\r = — 3 - + — 


V 2 y(r,0) 


d~\|/ 2 d\\f 1 


dr 


+ 


+ - 


r dr r 2 s i n Q d0 


3 ^sine^-) 

d0 


J 


or 


v = Vx 


lM) , 

rsinO 
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Orthogonality of \|/=const and ([)=const 


A contour on which \|/=const is called a streamline. A contour on which (])=const is called a iso- 
\potential line. It turns out that these two contours are 
orthogonal at any point in the fluid. To see this first note that 

v = V4> 

From the geometric meaning of gradient (see page 6 of 
Notes), we know that V(|) and hence v is normal to a 
(f>=const surface (see figure at right). Second, recall that v 
can also be written in terms of streamfunction as 


for R.C.C.S. 


v = V\|/xk 


Recall that the cross product is a vector which is orthogonal 
to two vectors being multiplied. Thus v, V\|/, and k are 
mutually orthogonal. Since v and V(f> point in the same 
direction, V(f> and V\| / must also be orthogonal. 



Streamlines, Pathlines and 
Streaklines 


Potential Flow Around Sphere. Lines lie 
in plane of page, which could be any 
<b =const. plane. 


Streamline - a contour in the fluid whose tangent is | 

everywhere parallel to v at a given instant of time. 

Path Line - trajectory swept out by a fluid element. 

Streak Line - a contour on which he ah fluid elements which earlier past through a given point in space 
(e.g. dye trace) 

For steady flows, these three definitions describe the same contour but, more generally, they are 
different. 
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Physical Meaning of Streamfunction 


The precise meaning of streamfunction is somewhat 
different for 2-D and axisymmetric flows. Let's focus on 
2-D flows normal to a cylinder (not necessarily with 
circular cross-section, axis corresponds to /-axis in 
R.C.C.S.). To motivate the somewhat lengthy analysis 
which follows, we first state the physical meaning. First 
we observe that material points follow trajectories which 
can be described as \|/=const. Three such trajectories are 
shown at right which he in the %y-plane. When these 
trajectories are (mathematically) translated along the z- 
axis a distance L they sweep-out \|/=const surfaces. 



No fluid crosses these surfaces: there are like the walls of a tube. Since no fluid leaves or enters this 
“tube”, conservation of mass means the mass flowrate must be a constant at any point along the tube. 
For an incompressible fluid, the volumentric flowrate is also constant. Suppose we denote the 
volumetric flowrate between any two of these \|/=const surfaces as AQ; then it turns out that 


A Q 

~Y = V2 “Vi 


Thus \|/ might be interpreted as the volumetric flowrate per unit length between this particular streaming 
surface and the one corresponding to 
\|/= 0 . 


Now let’s show this. Consider an 
arbitrary open contour (C) in the xy- 
plane, cutting across the flow. Next, 
consider the surface (A) formed by 
translating this contour a distance L 
parallel to the z-axis. The volumetric 
flowrate across A is: 



Q = j A n • \da 

where n is a unit normal to A. Since nothing changes with z, we choose a short segment of the contour, 
having length ds and of length L as our differential area element. 

da = L ds, 


The flowrate becomes: 


— = | n- v ds = | V\|/- dr 

c c 


(30) 
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where n is now normal to C and lies in the .i'y-planc. The key to this proof is deriving the second 
equality in the expression above. 

First, v can be expressed in terms of the gradient of the streamfunction using (27): 

v = V\|/xk (27) 

Both v and V\| / he in the plane of the page, whereas k is 
perpendicular to this plane (points out of page). Post -crossing V\|/ 
by the unit vector k does not change the magnitude but rotates V\|/ 
by 90° clockwise. If instead, we pre-crossed by k we would rotate 
V\|/ by 90° counter-clockwise. In either case, the cross product of 
k and V\|/ is a vector lying in the plane of the page and of the same 
magnitude as V\|/. 

The other term in the integrand of (30) 
is n ds, where d.s is the magnitude of a 
differential displacement along the 
contour, which we will call dr: 

ds = |<r/r| 

Since n is a unit vector n ds has the 

same magnitude as dr but is rotated by 90° . Both n ds and dr he in the 
plane of the page. Just as in (27), we can rotate one vector into the other 
by crossing with k: 

drxk=nds (31) 

into (30) yields 

n -v ds = v-n ds = (V\|/xk)-(<r/rxk) = V\\t-dr = d\\l (32) 

The 3 rd equality above says that the dot product of the two rotated vectors is the same as the dot 
product of the two vectors without rotation (since they are both rotated by the same amount). (32) into 
(30) yields: 



Substituting (27) and (31) 



Q 

L 


| n -\ds= | d\\i = \|/(g) -\|/(P) 

Cp-^Q Cp_>Q 
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To extract the physical meaning of this results, consider two 
contours, denoted as C 1 and C 2 in the figure at right. 
Notice that if C\ coincides with a streamline, the velocity is 
parallel to the contour at every point any no fluid crosses it: 

then: QIL = 0 

and \y(Q) = y(P) = y 2 

Thus \| / = const, along a streamline 


u - U| 



On the other hand, if the contour cuts across two 

streamlines (see contour C 2 in figure at right), then the difference in value of \| / corresponding to two 
different streamlines is just the volumetric flowrate of fluid held between the two streamlines (per unit 
length in the z-direction): 


A\) r = A QIL 


Incompressible Fluids 

By “incompressible fluid” we are usually referring to the assumption that the fluid’s density is not a 
significant function of time or of position. In other words, 

f + V.(pv) = 0 

can be replaced by V • v = 0 

For steady flows, dp/dt = 0 already and the main further requirement is that density gradients be 
negligible 


Identity C . 1 

V-(pv) = p(V-v) + v-Vp = p(V-v) 

Since flow causes the pressure to change, we might expect the fluid density to change — at least for 
gases. As we shall see shortly, gases as well as liquids can be treated as incompressible for some kinds 
of flow problems. Conversely, in other flow problems, neither gas nor liquid can be treated as 
incompressible. So what is the real criteria? 

For an ideal fluid (i.e. no viscous dissipation to cause V7), density variations come about primarily 
because of pressure variations. For an isentropic expansion, the compressibility of the fluid turns out to 
be: 
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1 

J2~ 


5 C 


where c = speed of sound in the fluid 

Thus changes in density caused by changes in pressure can be estimated as 


Ap ~ -r-kp (33) 

c 

According to Bemoulh’s equation, pressure changes for steady flow are related to velocity changes: 

2 

— + - — = const, or Ap = - Ap Av 2 (34) 

p 2 1 


(34) into (33): 


Ap 


p Av 
2 ,2 


Tlie largest change in density corresponds to the largest change in v 2 , which is v max 2 - 0 


( Ap 


_ i 

f \2 

v max | 

p 


~ 2 

c 

V K 

1 J max 

V J 


If the fraction change in density is small enough, then it can be neglected: 
Criteria 1 : v max « c 

for air at sea level: c = 342 m/s = 700 mph 

for distilled water at 25°C: c = 1500 m/s = 3400 mph 

For unsteady flows, a second criteria must be met: 


Criteria 2 : x » — 

c 

where x = time over which significant changes in v occur 

l = distance over which changes in v occur 
l/c = time for sound to propagate a distance / 
For steady flow x = °° and Criteria 2 is always satisfied. 
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Any fluid can be considered incompressible if both criteria are met. 

Viscous Fluids 

To resolve d’Alembert's paradox, we need to introduce another force into our force balance. This 
force can be thought of as resulting from friction between the fluid elements. Friction gives rise to 
viscous heating which represent an irreversible conversion of mechanical energy into heat. Indeed, 
friction is the main difference between an ideal fluid and a real fluid. 

friction — > irreversible deformation — > non ideal flow 

Tensorial Nature of Surface Forces 

Friction is a surface force like hydrostatic pressure, but unlike pressure, friction is not isotropic: 
isotropic - independent of orientation (direction) 

We say that pressure is isotropic because the magnitude of the pressure force is independent of the 
orientation of the surface on which it acts. Recall: 

dF„ = -n pda 

|t/FJ = pda 

which is independent of n (orientation). Viscous friction does not 
have this property: 

\dFJ{ depends on n 

but it's magnitude is proportional to da: 

I dFA °c da 

Thus it makes sense to talk about the force per unit area, which we usually call pressure. In the more 
general case in which the magnitude of the force per unit area might depend on the orientation of the 
surface, we use the name stress. In general, all surface forces can be lumped together and written as 

dF sur f = t(n) da 
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where t(n) denotes the surface force per unit area acting on the body through the surface element whose 
orientation is given by the unit outward normal n. We call this quantity the stress vector.* In what 
follows, I will try to convince you that the effect of orientation of the surface can be expressed as 

t(n) = n • T 

where I is called the stress tensor. While I depends on position and possibly time, it does not depend 
on the orientation of the differential surface element da, which is given by n. T is sometimes call the 
“state of stress” of the fluid. 

First, let’s try to understand better what we mean when we say a 
material experiences some “stress”. Consider a block of some material 
which bears some externally applied equilibrium load (see figure at 
right). By “equilibrium” I mean there is no net force and no net torque 
applied to the body: thus it is at mechanical equilibrium and has no 
tendency to accelerate. The material might be a fluid or a solid, but what 
we are about to say is easier to imagine if we think of material as a 
solid block. 

Now consider some mathematical surface inside the material 
(indicated by the dotted line in the figure above). What are the 
forces exerted across this mathematical surface? To answer this 
question, suppose we actually separated the block into two pieces 
by physically cutting along this surface without changing the loading. 

Block “1” now experiences an unbalanced load -F while block “2” 
experiences an unbalanced load +F. Once separated, both blocks 
would tend to accelerate in opposite directions. 

Why don’t the two halves accelerate when connected? Apparently, half “2” must have exerted a force 
on half “1” which we denote as t(iq)A, where A is the area of the cut face, and which equals 

t(n | )/V = +F 

while half “1” exerted an equal but opposite force on half “2” 

t(n 2 )A = -F = -tfn^A (35) 

When expressed per unit area, this internal force between the two halves (when they are physically 
connected) is what we mean by the “stress” experienced by the material under load. Given that n 1 = 
-n 2 , (35) tells us that 


*see Whitaker, Introduction to Fluid Mechanics, Chapt. 4. 
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t( n 2 ) = -t(-n 2 ) 


More generally, for any mathematically surface having orientation n, we can write t(n) = -t(-n) 

or t(-n) = -t(+n) (36) 

In essence, this is just a statement of Newton’s Third Law (for every 
action, there is an equal but opposite reaction). 

Now let’s generalize to three dimensions. Suppose we know the 
distribution of stress (i.e. the “loading”) on all six faces of a block of 
material. Furthermore, suppose this loading is an equilibrium loading 
(no net force and no net torque on the block). Let’s try to calculate 
the stress on a mathematical surface cutting through the material at an 
arbitrary angle. Let the orientation of the mathematical surface be 
given by the unit normal n. 

Problem given the surface stress on mutually perpendicular planes 
[i.e. given t(i) for v=const plane, t(j) for y=const plane, and t(k) for 
z=const plane], calculate the surface stress on a plane of arbitrary 
orientation specified by the unit normal n. 

Given: t(i), t(j), and t(k) 

Find: t(n) 

Solution: We choose the tetrahedron ABCO as our “system.” For the surface forces to be balanced* 

|t(n)<r/a = 0 
A 

We evaluate this surface integral by subdividing the surface A into the four faces of the tetrahedron: 
planes ABC , BCO. AOC, and ABO. For each surface, we need to evaluate the outward normal n, the 
stress vector t and the surface area A. The results are tabulated below: 



*More generally, in our force balance, we should include body forces and inertia as well as surfaces 
forces. If we let V— >0 then the volume integrals for body forces and inertia vanish more rapidly than 
surface forces. Thus for the surface A enclosing a tiny volume V, we must still require ^ t da = 0 for the 

A 


surface of every differential volume element. 
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Plane 

Outward 

Normal 

Area 

Stress 

Vector 

ABC 

n 

A 

t(n) 

BCO 

-i 

n- iA 

t(-i) 

AOC 

j 

njA 

t(-j) 

ABO 

-k 

n • kA 

t(-k) 


Applying the Mean Value Theorem to our surface integral: 

j t da = 0 = [(t(n)) + (n • i)(t(-i)) + (n • j)(t(-j)) + (n -k)(t(-k))] A 
A 

Dividing by A cancels the A out of our expression for the integral. Taking the lim it as all the dimensions 
of the block vanish (i.e. as A— >0) allows us to replace the unknown averages with their limit, which is 
the value of the vector at the point that the tetrahedron collapses about. 

t(n) + (n • i)t(-i) + (n • j)t(-j) + (n- k)t(-k) = 0 

or t(n) = (n • i)t(-i) - (n • j)t(-j) - (n • k)t(-k) (37) 

Next we apply (36). Of course, we can replace n in (36) with i or j or k. Using (36) in (37): 

t(n) = (n- i)t(i) + (n • j)t(j) + (n- k)t(k) 

Recalling the definition of dyadic product, we could re-write this expression as 

t(n) = n • it(i) + n • jt(j) + n • kt(k) 
t(n) = n • [it(i) + jt(j) + kt(k)] = n ■ T 

The sum (inside square brackets) of these three dyadic products is some second-rank tensor, which is 
independent of n. We denote this tensor by T. 

Significance : to calculate the surface force on some differential surface area da having orientation n, we 
just multiply this expression for the stress by the area: 

d¥ S urf =n-Tda (38) 

where I is called the stress tensor and dF sur f represents the net surface force (all contributions) acting 
on the body whose outward normal is n. While T does not depend on n, it might depend on position 
inside a solid which is nonuniformly stressed; thus 

T = T(r) 
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represents the state of stress of the fluid. In terms of the components of the tensor, we say that 

Tjj = / lh component of the force acting on a r ( - = const surface 

To make this more understandable, it might help to express I in terms of other variables in a familiar 
problem, say an ideal fluid. 

For an ideal fluid: J = -pi 

To show that this is correct, we will compute the surface force and show that it reduces to the fa mil iar 
force due to pressure: 


d¥ sut f = n • T 'cla = n • (-pl)da = -pin • I )da = -pnda = dFp 


Generalization of Euler's Equation 

Recall that Euler’s equation was derived by applying Newton's Second Law ( F=ma ) to any fluid 
element. To generalize this result to include friction, we replace the “pressure force” by the surface 


force. 


Instead of: 

r Dy 

1 dV=¥g+Fp 
V 

we have: 

r Dv 

J P p )l ~ Fg + ^ surf 


V 


V 


V 


Applying the Divergence Theorem and combining the three volume integrals: 

)dV = 0 


r f Dy ^ 

jl P — -Pg-V-T 


V 


Dt 


Since this must result for all choices of V in some region, the integrand must vanish at every point in that 
region, or: 


Dy 1 

— = g+-VT (39) 

Dt p = 

which is a more generalized version of Euler's Equation. 

For an ideal fluid: J = -pi 
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identity 

C.ll 

and V • T = - V • ( pi ) = -Vp 


For a proof of identity C.ll, see Hwk #3, Prob. 4. For a real fluid, there is an additional contribution 
to T from friction: 

For a real fluid: T = -pi + x (40) 

where x is called the viscous stress or the deviatoric stress (since it represents the deviation from 
ideal). Note that: 

CAUTION: x in these notes (and in Whitaker) — > -x in BS&L 

By convention, Xp in these notes is positive for tensile stresses (which result from stretching a solid rod) 
whereas Xp in BS&L is positive for compressive stresses (which result from putting a fluid under 
pressure). Although BS&L's notation might make more sense for fluids (which usually do not 
experience tensile stresses), we will use the other convention because it is more commonly used in 
continuum mechanics: in particular, this is the convention used by Whitaker. 

When written in terms of scalar components, the above equation between T and x represents 9 
equations in 10 unknowns (the nine components of x plus p). Clearly p and x cannot yet be uniquely 
determined, given the state of stress T. One additional relationship is required. For a real fluid, we 
somewhat arbitrarily define p as the average of the three diagonal components of T: 

p = 4(T:l) (41) 


where 


T: I = 


i j 




or, in Cartesian coords.: T : I = T xx +T yy +T zz 

In any coordinate system, T : I is called the trace of T. Then p can be thought of as the isotropic 
contribution to stress (that part of the normal stress which acts equally in all directions) while x 
represents the remainder, or the nonisotropic part. (40) might be regarded as decomposing x into an 
isotropic part and a nonisotropic part, which is a common thing to do with tensors. 


The choice of p made in (41) also represents the thermodynamic pressure ; i.e. the hydrodynamic 
p is now the same as the p appearing in the thermodynamic equation of state: 


p = p(p,D 


To summarize, we decompose the state of stress into two contributions: an isotropic pressure: 
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and a deviatoric stress: x = I - (1/3)(X : 1)1 

Taking the divergence of (40): V • T = -V/; + V • t 

and substituting it into Euler's Eq. (39): 



Momentum Flux 

Up to this point we have spoken of I as the stress on the fluid. I can also be thought of as 
momentum flux. To convince yourself that momentum is being transported (like heat and mass), 
consider the problem of unsteady simple shear flow. At time 1=0. an initially quiescent fluid confined 
between i nfi nite parallel plates is disturbed by imposing motion on the upper plate. The velocity profile 
gradually develops into linear shear flow. 

Note that, like all other fluid elements, the fluid ► U 

element at y=hl 2 undergoes acceleration: 

dvjdt > 0 at y=hl 2 

In other words, the fluid element is gaining 
momentum. How can it acquire momentum? 

Answer: momentum is transported to y=h/2 from 
above through friction between fluid elements. There 
are two directions associated with transport of 
momentum: 


1) direction of the momentum being transported 
(in this example, x -momentum is transported) 

2) direction in which the momentum is transported 

(y) 

For this reason momentum flux must be a 2nd rank 
tensor. It turns out that: 
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-X = diffusive flux of momentum* 
pvv = convective flux of momentum 

In multicomponent systems, the flux of any species i due to convection is written as the product of 
concentration and velocity: 

conv. flux of species i = c{v [=] mol-cnr—s -1 [=] rate/area 
where c t [=] mol/cm 3 

More generally, by “flux,” I mean: 

flux — that tensor (of whatever rank), which when pre-dotted by n da, gives the rate of 
transport through the surface element having area da in the direction of its unit normal 

n. 

For example, the flux of fluid mass by convection is 

pv [=] g-cnr 2 -s _1 [=] rate/area 
Proof : pre-dotting by n da we obtain: 

(n da) • (pv) = pn • \da = p (dq) [=] g/s 

which represents the mass flowrate across da. Indeed the convective flux of total fluid mass also has the 
form of concentration times velocity since 

p = concentration of mass [=] mass/vol. 

It might come as less of a surprize that pvv is the convective flux of momentum if you rea liz e that 

pv = cone, of momentum [=] momentum/vol. 

After all, momentum is mass times velocity. So 

(pv)v = pvv = conv. flux of momentum 

Example: Apply conservation of linear momentum to an arbitrary fluid system. Thus prove that T and 
pvv are momentum fluxes, as claimed. 

Solution: We choose a system that has fixed boundaries in the laboratory reference frame: in other 
words, V ± V(l). Referring to the discussion above in which pv is the concentration of momentum, pvv 


* "Diffusive" means it results from random co lli sions between molecules. 
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is the momentum flux by convection, and -I is the momentum flux by diffusion, conservation of 
momentum requires 


— J p \ clV = -jin-pvvdtf +Jpg clV 

14 , 2 43 2 4 3 1^2 4® 142 <© 

accumulation in by convection in by diffusion in by external forces 

Applying the divergence theorem combining the resulting volume integrals, and invoking the result for 
arbitrary V, we would obtain 

-VT-pg=0 

[v-(pv)Jv+pv-Vv 

Expanding the partial time-derivative of the product and expanding the divergence of the dyadic product 
using identity C.8 


U p ; ] 

Bp Bv 
— v+p — 

a t dt 


+ 


^+ v -(pv) 

L r442 4 43 

o 

Finally, we recognize the factor inside square brackets must vanish according to the general continuity 
equation. After this term is dropped, the remainder is Euler's equation: 

dv 

p — + pv-Vv = V-T + pg 
at — 

Comment: In writing the statement of conservation of momentum, we have a term representing 
transport of momentum into the system by the action of external forces. Clearly, an object falling from 
rest in a gravitational field acquires momentum through the action of gravity. Does this acquisition 
represent transport to the body from outside, or does it represent a "generation" term? If it is 
spontaneous momentum generation, then momentum is not conserved. 

With a little reflection, we can convince ourselves that the action of gravity is “transport” and not 
“generation.” Some of the earth's momentum is being transported to the fading object. When the object 
eventually collides with the earth and comes to rest again, that momentum is transported back to the 
earth. The total momentum of the universe has not changed: momentum is conserved. 

Actually, the term we call “diffusion of momentum” also arises from the action of an external force: this 
time it's the action of “surface” forces, rather than “body” forces. 


dv 

v +p — + pv-Vv- V-T-pg = 0 
at ~ 
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Response of Elastic Solids to Normal (Uniaxial) Stress 

By generalizing, the force balance plus the mass balance now contain more unknowns than 
equations: 


x, v, and p — > 9+3+1 = 13 unknowns 
Euler + Continuity -+3+1=4 equations 

Missing is the constitutive equation which is an empirical description of how the fluid or solid material 
responds to stress. Obtaining this relationship is an important objective of that field of science known as 
continuum mechanics. 


Let's start by considering solids, whose response is more familiar. Suppose I try to stretch a rectangular 
bar by applying a tensile force, F. 

We will attempt to describe the response of 
the bar under conditions of mechanical 
equilibrium. To have the bar - stretch instead of 
accelerating as a result of the force, I must 
apply an equal but opposite force to the other 
end. Let the v-axis be aligned with the 
direction of the applied force: 

F x = |F| 

At equilibrium, the length of the bar will increase by an amount 5 X . Hooke's law tell us that: 

F I 

5 

II 

L ^y L -‘z 



but 


L y L z 


is the applied stress. The two subscripts on stress denote the two directions associated with it: the first 
subscript denotes the orientation of the surface the force is applied to (v=const) while the second 
denotes the direction of the applied force (x). 

Since the deformation is proportional to the original length of the bar, it makes sense to define the 
deformation per unit length, which is called: 
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strain : 


'XX 




Later we will see that strain is a second-rank tensor. For now, we will interpret the order of the 
subscripts as follows: the first subscript gives surface (in the above case, an v=const surface), while the 
second subscript gives the direction of the deformation (in the above 
case, the x- -direction). 

Experiments show that the strain increases with the applied stress as 
shown in the figure at right. At low stress levels, the strain is directly 
proportional to the applied stress: Hooke's law* for purely elastic 
solids is: 

T = F? 

i XX ^^xx 

E [=] force/area 

where E is called the Young's Modulus. It turns out that deformation also occurs in the y and z 
directions: 



e yy ~ £ zz ~ ~ ve xx 

where v is called Poisson 's ratio. For most materials 

0.28 <v <0.33 


Response of Elastic Solids to Shear Stress 


Instead of a normal force, suppose I apply a shear force on the upper face. To keep the object 
from translating, I must apply an equal but opposite force on the lower 
face. 


This generates a force couple or torque, which will cause the body to , r 
spin. To prevent a steady increase in rotation speed, I must apply an - 1 
equal but opposite torque. Recalling that torque is force times lever 
arm: 


„ ‘ 


\ 

'± l 3 




FI -FI 
r x^y ~ r y^x 


* Robert Hooke (1635-1703). English experimental physicist. Hooke’s law first stated in 1660. 
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Dividing by L x L y L z \ 


Ex 

L X L Z 


LyL 


z 


or 


T =T 

*yx ± xy 


where T yx = .r-component of force/area acting on y=const surface. Note that this imphes that the stress 
tensor is symmetric. This turns out to be true for virtually all loadings.* 

This loading produces a shear deformation in which 

Ex “ L y 


and 


5 v “ L x 


Hie corresponding definition of shear strain is: 

f o s y 


8 8 

.1 

\ L y L x j 


— £ = £ 

° xy ° yx 



Since we applied a symmetric stress, by applying forces in both the x and y directions, we average the 
strains in the x and y directions to obtain a symmetric strain. Moreover, this definition yields a strain 
which is invariant to rotation of the xy axes. As with normal stress, shear stress produces a shear strain 
in direct proportion to the stress: 


T xy = 2 Ge xy (pure shear) 

where G is call the modulus of elasticity for pure strain. Although the units are the same, the value 
of 2 G is different from that of E. 


Generalized Hooke's Law 


Now let's try to generalize to some arbitrary loading which might involve both shear and normal 
stresses. 


Consider two arbitrary material points in 
the material. Let x denote the position of 
the second material point relative to the 
first, before the load is applied. After the 
load is applied, both material points 


rc-Jalivc 
displacement 
8— of two material 
points by load 



* For a proof based on the assumption of local equilibrium, see W:4.3 
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move as shown by the dashed lines in the figure at right. After the load the position of the second 
material point relative to the first becomes x+8 , where 

8 (x) = relative displacement 

of two material points initially located at x. The resulting strain can be generalized as: 


or 

Note that this general definition 
normal or pure shear stress: 




V8 + (V8) f 


2 


V dxj a Xj J 


of strain reduces to earlier expression for strain in the cases of pure 


unixial normal stress: 


£ = — 

Yxx 2 


as 


X + 55 x ^ 


a.v c)x 


s Y 

= -*- and T 
L v 


- F F 

XX ^ c xx 


pure shear: 


'xy 


1 

P 8 .v , 38 ,j 

_ 1 

A +57 

2 

dx dy , 

2 

v L x L y y 


and T xy = 2Ge xy 


If the strain components are all small, we might reasonably suppose that Hooke's law generalizes into a 
linear relationship between any component of strain and the nine components of stress (or vice versa): 


3 3 

% = EE Cijkl ^ kl 

k = 1 / = ! 


(42) 


There are then nine coefficients for each component of stress, making a total of 81 possible coefficients. 
But just by requiring: 


• symmetry of I (T[-=T- n ) 

• isotropy of material (e.g. same Young's modulus applies to uniaxial stress in x -direction as 
for y and z.) 


it can be shown that the number of independent constants is reduced from 81 to only 2. It is customary 
to express the stress tensor as the sum of two tensors: one isotropic [(§ : j)J] and the remainder [g - 
(g : 1)1]. Denote the two independent constants as k \ and k 2 and use them as weighting factors in the 
sum 
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T = *^-(|^]+* 2 (g^ 

remainder isotropic 

This can be rearranged T = ky £ + (k 2 - ky : ijl 

Sometimes the constants are choosen as the two coefficients in this new form 

T = 2r|e + ?i(£:l)l 

which represents the generalization of Hooke's law of elasticity. Thus of G, E and v, only two are 
linearly independent properties of the material; using the above relationship we can show (HWK Set #5, 
Prob. 1) that 


v = 



Response of a Viscous Fluid to Pure Shear 


Suppose I have a thin layer of fluid held between parallel plates and I apply a force F x to the upper 


plate. To cause deformation, rather than simple 
translation of the system, I must apply an equal but 
opposite force to the lower plate. 

As before, this produces a force couple or torque. 
To keep the system from rotating, I must apply an 
opposing torque. Once again this equilibrium state 
corresponds to a symmetric stress tensor. From our 
experience, we know that eventually the upper plate 
will slide past the lower plate at some steady relative 
speed U x . 



F / 

1 x 

response: U x °c 

L x L z 


This speed represents a rate of deformation: 


U x 


dd x 

dt 


and the speed per unit thickness represents a rate of strain : 
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<r/8 


Ux _ 

/dt 

_ d 

Ly 

h 

dt 




y L yj 


cle 


yx 


dt 


Finally, the applied force divided by the area over which it's applied represents a shear stress: 

F . 

— — = 7 = T = T 

, , ± yx L yx L xy 

L X L z 


Rewriting the proportionality as an equality, we have: 


de 


^ yx “ M-' 


yx 


dt 


which is called Newton 's Law of Viscosity (1686). Alternatively, BSL point out that 

dv Y 


U r dv Y 

~ — — and 


dy 


( x -^)bsl ( Xj ?) 


= -M-- 

notes dy 


Generalized Newton's Law of Viscosity 


This generalization of Newton's law of viscosity to arbitrary loading parallels that of Hooke's Law 
with the strain tensor replaced by the rate of strain. The main difference is that now the deformation is a 
function of time: 


8 (xj) = displacement of material pt. 

= deformation 

The trajectory of a material point initially located at 
x(0) is given by: 

x[x(0),fj = x(0) + 8 [x(0),t] 


before load 



(lisplartfinenl 
SD) - of material point 
by load 


after load is 
applied for time t 


Keeping the material point fixed (i.e. x(0) is constant) 

while we take the time derivative is the same as taking the material derivative of position. Thus the rate 
of deformation of fluid elements is: 


J8 _ Dx _ 
dt Dt 


(rate of deformation) 


given by the fluid velocity. Just like we define the gradient of the deformation to be the strain: 


e = i[v8 + (V8) r ], 


(strain) 
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the gradient of the rate of deformation must be the rate of strain : 


d = 


Z)£ 

Dt 



j V v + (Vv) / j (rate of strain) 


Newton’s law of viscosity is a linear relationship between the stress and the rate of strain. Generalizing 
that linear relationship again leads to a relationship like (42), except that the strain tensor £^/ is replaced 
by the rate of strain tensor d^. Once again, in order for the stress tensor T or x to be symmetric and 
for the material to be isotropic, only two of the 81 coefficients Ck-^/ are independent. We use the two 
independent coefficients to multiply the isotropic part of x and its remainder: 

x = 2pd + (k - -j |_i)(d: l)l 


where p is the usual viscosity and K is called the second coefficient of viscosity, d : I has a special 
significance, which we will now point out. As we showed on page 70, d : I is the trace of d: 

d : I = dy \ + 6?22 ^33 

Writing in terms of the velocity v: d[j=\ 


OV; °v 


-+■ 


J 


d.X ; dX; 


SO 


d 11 


dv 1 


- + 


dx^ 


dvi 


_■ T dv! 
d: I = L 


+ - 


dv-> dv 


3 _ 


— d.i'[ d.i '2 dxi 


= V- v 


So the trace of the rate of deformation tensor is just the divergence of the velocity. Newton 's law of 
viscosity becomes: 


x = 2pd + (K-|p)(V-v)l 


For an incompressible fluid: 


V-v = 0 


leaving 


x = 2pd = p 


Vv+ (Vv) ? 


(incompressible) 


Navier-Stokes Equation 

Once again, let’s go back to Euler's equation, generalized to account for the tensorial nature of 
viscous friction: 
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pD\/Dt = pg - Vp + V • i 
For an incompressible fluid, V • v=0 and 

I = 2pd = p[Vv+(Vv) f ] 

V • x = p[V • Vv + V • (Vv) ? ] 

but using identity C.10: V • (Vv) f = V(V • v) = 0 

leaving V ■ x = pV • Vv = pV 2 v 

D\ 2 

Euler's equation becomes: p = pg - \p + pV~ v (p ,p=const) 

which is known as the Navier-Stokes Equation (1822).* Now we have as many equations as 
unknowns: 


v, p — ^ 4 unknowns 
N-SE, Continuity — > 4 eqns 


Boundary Conditions 


But to successfully model a flow problem, we need more than 
a sufficient number of differential equations. We also need 
boundary conditions. 

A typical boundary is the interface between two immiscible 
phases — either two fluids or a fluid and a solid. One such 
boundary condition which can generally be applied is the no slip 
condition : 



y i= y ii 


* Sir George Gabriel Stokes (1819-1903): British (Irish bom) mathematician and physicist, known for 
his study of hydrodynamics. Lucasian professor of mathematics at Cambridge University 1849-1903 
(longest-serving Lucasian professor); president of Royal Society (1885-1890). 
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For example, in the problem of uniform flow around a stationary solid 
sphere, this requires: 

at r=a : Vy= v v = 0 

which means that there is no flow of fluid across the boundary 

at r=a : v r = 0 

We also assumed this boundary condition when we solved the potential 
flow problem. But “no slip” also means: 

at r=a : v e = 0 

Note that our potential flow solution did not satisfy this second equation: 

for pot’l flow: v r (a,0) = 0 

v e (<3,0) = -(3/2)£/sin0 * 0 

In addition to d'Alembert's paradox, potential flow fails to satisfy the no-slip condition.* 

For a fluid-solid interface, in which the velocity of the solid phase is known, the no-slip condition is 
sufficient. But in the case of fluid-fluid interface, the velocity of the second fluid is usually unknown. 
Then no-slip just relates one unknown to another. 

A second boundary condition can be obtained by considering the stresses 
acting on the material on either side of the interface. Suppose we were 
to apply a loading as on page 66 to a two-phase region straddling the 
interface, as suggested in the figure at right. Note that the loading is 
balanced: that is there is no net force on the system. 

If we were to split the system into two parts along the 
interface, each of the two halves would tend to accelerate. 

This suggests that, when the two phases are in contact, each 
exerts an “internal” force on the other, as shown in the figure 
at right. These forces are equal but opposite: 

^lon2 — " ^2onl — ^ 


* Although "no slip" is usually applicable, there are at least two situations where no slip might not be 
applicable: 1) when the mean-free path of gas molecules is comparable to the geometric dimension, and 
2) when a liquid does not wet the solid. 
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which might be thought of as “Newton’s Third Law”: for every action, there is an equal but opposite 
reaction. Using (38) (on page 68) to express the forces in terms of the stress tensor: 

112 ■ Tj A = — n | -T A = F 

where A is the area of the interface and where n, is the normal to the interface pointing out of phase i. 
However, from the geometry, we can deduce that 

n t = -n 2 = n 

Newton’s Third Law becomes: n X| = n - T 0 (43) 

When the interface is highly curved (e.g. a small oil droplet in water), then surface tension can produce a 
discontinuity in the normal components of the above force, which has not been included in the above 
analysis [see L&L, Chapt. 7 or Hunter, Vol I., p237f]. The more general form of this boundary 
conditon is 


n ‘CTi - T 2 ) = y(V 5 • n)n 


(44) 


where y is another property of the fluid called the surface tension and 

V lV = (I - nn) ' V 

is the surface component of the V operator. We will have more to say about this near 
the end of the course. For now, we will neglect surface tension effects. 

V v • n = curvature of surface [=] m _1 

For a flat surface, n is independent of position along the surface so that V v • n = 0 and 
(44) reduces to (43). 


Exact Solutions of N-S Equations 

Exact solution of the Navier-Stokes equations presents a formidible mathematical problem. By 
“exact” I mean: 

exact — neither viscous nor inertial terms are neglected (i.e. approximated by zero, as opposed 
to being identically zero) 

One difficulty is the non-linear inertial term. Most of the powerful mathematical techniques (such as 
eigenfunction expansions, used in “separation of variables”) only work when the P.D.E. is linear. 
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Of course, if we can neglect the higher-order viscous terms, then we can cope with the non-linearity 
using a velocity potential, as we did earlier in solving problems in potential flow. However, viscous 
terms are seldom completely negligible and leaving them in the equation makes the problem much more 
difficult by increasing the order of the P.D.E. Nevertheless, it is possible to find exact solutions in 
certain cases, usually when the inertial terms vanish in some natural way. We will now examine a few of 
these problems having exact solutions. 


Problems with Zero Inertia 

First, let’s consider problems in which the fluids elements travel along straight streamlines at constant 
velocity. Then their acceleration vanishes identically. 


Flow in Long Straight Conduit of Uniform Cross Section 

Suppose we have pressure-driven flow in a long straight 
conduit whose cross section does not vary along the flow. In 
mathematical terms, the conduit is a cylinder of arbitrary cross- 
section. Define Cartesian coordinates such that the axis of the 
cylinder corresponds to the z-axis. In a very long pipe, we expect 
that v z will depend on z (as well as x and y) near the entrance and 
exit of the pipe. 



\ 


M ►+* ►H M 

Entrance Fully Exit 

region developed region 


In particular, near the entrance, we say the velocity profile is “developing”; i.e. evolving with z. Outside 
the entrance and exit regions v, will be independent of z. This situation is called fully developed flow. 
For fully developed, steady flow: 


v z = v . (x >y) 
v x = v y = 0 

Note that this automatically satisfies the Continuity Equation: 


V • v = dvfdz = 0 
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For this form of solution, it turns out that the nonlinear inertial terms automatically vanish: 


v Vv = 0 


To convince yourself of this, consider steady flow. Since the velocity of a fluid element is constant along 
a straight streamline: 

for steady flow: v • Vv = Dv/Dt = a = 0 

In other words, fluids elements are not accelerating. Thus inertial forces are identically zero in steady 
pipe flows. Strictly speaking, Reynolds number should not be thought of the ratio of inertial to viscous 
forces in this problem, since inertia is zero for laminar flow although the Reynolds number is not zero. 

For a steady flow, the N-S equations are: 

0 = pV 2 v - Vp + pg 

Note that for a vertical pipe:* p = p(z) 

but more generally for a horizontal or inclined pipe: 

P =p(x,y,z ) 

where the dependence on position in the cross section arises from the contribution to pressure from the 
hydrostatic head (i.e. from g). For this and some other problems, it’s helpful to decompose the total 
pressure into contributions from gravity (i.e. hydrostatic pressure, p h ) and from flow (called the 
dynamic pressure , P ) 


p = p h + P 
Vp = Wp h + VP 
pg 

From our earlier analysis of hydrostatic equilibrium (see page 39), we know that Vp h = pg . 

Note that Vp - pg = VP 

Next, we substitute this into the N-S Equation. Expanding them in component form: 


* In the absence of gravity, dp/dx and dp/dy must vanish: see the x- and y-components of (45). 
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- component: 
y - component: 


0 = 
0 = 


dP_ 

c)x 

dP_ 

dy 


\p = P(z) 


z-component: 


0 = (iV 2 v, - — 

z dz 


(45) 


where we have written the last term as the total derivative instead of the partial derivative because the 
first two equations require that P be a function of z alone. We can immediately deduce that P must be a 
linear function of z 


^ 2 $ 

g{x,y ) 



since the dP/dz is independent of x and y, while the velocity profile is independent of z, so (45) requires 
that the two functions of different variables be equal: 


g(x,y) =/(z) = const, w.r.t. x,y,z 


which can only be hue if both functions equal the same constant: thus P(z) must be a 
linear function. For steady flow in a circular conduit of radius a, the “no-slip” b.c. 
requires 

at r = a: v z = 0 (46) 



Since neither the b.c. nor the differential equation contain any dependence on 0, we 
expect the solution to be axisymmetric about the z-axis: 


(45) becomes: 


v z = v z ( r ) 



r dr v 


dv z ' 
dr j 


1 dP 

= const 

(l dz 


The general solution of this equation is: 


v 


z 


(r) 


1 dP 2 

r + ci In r + C 2 

4(1 dz 


Requiring the solution to be bounded at the center of the tube (as r— >0) forces us to choose c 1 =0 while 
the remaining constant can be chosen to satisfy the no-slip condition (46). The particular solution is the 
familiar parabolic velocity profile: 
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/ \ 1 dP ( 2 2\ 

v 7 (r) = \r -a 

ZW 4[Ldz { 1 

The volumetric flowrate through the conduit 
is computed from: 

Q = jn-\da = J a v,(r)(27tr)<ir 






A 


or 


Q 


7 la 


8|i 


side view 

4 r 

dz 



da 

(area of 
shaded 
region) 


end view 


which is called Poisueille's Formula. This formula was been derived for a number of different cross 
sections. In general 


Q 


kA A ( dP 3 




dz ) 


where 


A = cross-sectional area of duct 


and where k is some constant which depends on the shape of the duct; e.g. 
circle: k = l/87t = 0.0398 


square: 

ellipse: 


k = 0.0351 
8 

47t(l + £ 2 ) 


where e = b/a < 1 is the ratio of the minor to major axis. 


Flow of Thin Film Down Inclined Plane 

Suppose we have fluid overflowing some reservoir and 
down an inclined plane surface. Although there might be 
some entrance or exit effects (at the upstream and 
downstream ends of the plane), if the plane is sufficiently 
long compared to these regions, then what we see in 
experiments is a region in which the fluid flows downward as 
a film of uniform thickness. Let's try to analyze this central 
region in which the film in uniform. 
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Let the .r-axis be oriented parallel to the inclined plane in the direction of flow and let the y-axis be 
perpendicular to the plane. One of the boundary conditions will no “no-slip” at the solid surface: 

at v=0 : v x = 0 

At the very least, there must be flow in the x -direction: otherwise there could be no viscous force to 
balance gravity. In addition, the x -component must vary with y : v x must vanish at y=0 (no slip) and be 
nonzero for y>0. The simplest form of solution which is consistent with those constraints is: 

v = v x (y)e x 

Note that this form automatically satisfies the requirement of continuity for an incompressible fluid: 

V • v = dvjdx = 0 

Substituting this velocity field into the Navier-Stokes equations: 

x : 0 = -dp/dx + [ufiv jdy 2 + pg A . (47) 

v: 0 = -dp/dy + pg y (48) 

z: 0 = -dp/dz 

where g x = g • e r = gcosp 

and g y = g • e v = -gsinp 

Integrating (48) with respect to y: 


p(x,y ) = - pgysinp + c(x) (49) 

where the integration constant might depend on x, but cannot depend on z (according to z component 
above) or y. Now let's turn to the boundary conditions. Continuity of stress across the interface 
yields: 

at y-8 : n ■ T liq — n ■ n ' (~Pgasl =gas - ) 

Now the viscosity of air is about 0.001 times that of water. Then it is reasonable to treat the air as an 
inviscid fluid: i.e. neglect l gas . This leaves 


n '=liq ® Pgas 11 P atm 

where p atm = latm. Now e v is the unit normal to the interface in this problem 

e y * lliq = e y ’ ('Pi + l) = 'P e y + ^yx e x + %y e y + hz e z 


( 50 ) 
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When v x = v x (y), the only nonzero components of the deviatoric stress tensor are x fv =x vr * Dropping 
the other terms, (50) becomes: 


-p£y + tyx^x — ® yPatm 

Equating separate components: 

* -component: x yx (x,y= 5) = 0 (51) 

y-component: p(x,y=8 ) = p atm 

With this result, we can evaluate the integration constant in (49): 

P(x,y=8) = -pgSsinp + c(x) = p atm 
Thus c(x)=p atm + pgSsinp 

(49) becomes: p(x,y) = p atm + pg(5-y)sin|3 

Thus dp/dx = 0 and (47) becomes: 

pcl 2 v x /dy 2 = -pgcosp 

No-shp at the wall requires: 


at y=0 : 


v x = 0 


whereas if the stress in (51) is evaluated using Newton's 
law of viscosity, we also require: 


at y=5 : x yx = p dvjdy = 0 

Using these two boundary conditions, the velocity profile 
can be uniquely determined: 


v x{y) 


pgS cos(3 


2p 


(A 

UJ 


vS j 



For a Newtonian fluid, the 9 scalar components of the stress tensor are expressed in terms of the 
derivatives of the velocity field on pi 45 of Whitaker. These expressions (except for a change in sign) 
can be found on p88 of BSL or at the website for our course — 

http ://w ww . andre w .emu .edu/course/06-7 03/NLV_RCCS .pdf. 
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By integrating over a plane perpendicular to the flow, we can evaluate the flowrate: 

W 5 

2 = Jn-vr/n= J J v x (y)dydz M ^ 

A 0 0 I 


pglTS 3 cos [3 


where W is the width of the plane. 

end view 

Time Out : The main novelty of this problem is the 
treatment of the free surface. We treated the air as if it 

was inviscid, although it has some viscosity. How important is the drag imposed by the air? 
subject of Hwk #6, Prob. 2 


To answer this question, let’s consider a vertical fi lm of water in 
contact with a vertical fil m of air, as shown in the sketch at right. 
Let’s re-solve the problem and see how large 8 a must be for a given 
S H , before we can neglect the effect of the air. For fully developed 
flow, the velocity and pressure profiles should have the form: 


Ft =Lc(y) 

V v = V, = 0 


p(x,y,z) = 1 atm 



NSE X becomes: 


for water: 


„ _ d 2 v w x 

M ‘W T P wS 

dy 2 


for air: 


0 = (l 


a o 

dy 


where we have taken p a = 0. Applying “no slip” at each of the three interfaces 


at y=0: 


v x (o) = o 




at y- 5 W ' + 5 a : 


( 8 W + 8 fl )-0 
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The particular solutions to the two ODE’s are 


v]TM = |^(s w -y)y + ^ 

8 v 


:w=^ 


8 w + 8 « - y 


Hie interfacial speed U is unknown, but is choosen to match the shear stress at the interface: 


aty=5 • 


<x{K)=<x{^w) 


For Newtonian fluids, the stresses can be related to velocity profile: 


Using the velocity profiles of (52) and (53), this stress matching yields: 


is U U 

M'a c- 


Solving for U : 


U = ± 


t pA 


2 Evv , V-a 

8 w 8 « 


The flowrate of the water is 


Q _ fw/\ i _ Pw’£ 8 w , 1 t — Pw& 8 w 8 


/| M'w | Pa 


-= \<(y)dy = 


’jvS^w , 1 _ vv 

12(1 w 2 1 2(_i w M-w | M-q 

8 w 8 « 


If the air film is very thick, the flowrate becomes 


— ^ lim \Q}=^^ 
W 8 a ^ l ^ S 3(1 w 


which is the same expression we had in the Notes when the viscosity of air was completely ignored. 
Dividing (54) by (55): 
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/j M-w | M-a /j M-w | 

G = 1 5 w 5 fl = 1 8 « 

Goo 4 Fh L+ (^ l 4 Fh L+ 5^ 

8 a M a 8 fl 

When the viscosity of water is 1000 times larger than air (i.e. p w = 1000 p fl ), this gives 

4000 + ^ 

_G_ = i <k_ 

4 1000 + ^ 

8« 

To reduce the flowrate by 1% means = 0.99 , for which the air film thickness must be 

Goo 

S fl =0.074 8 w 

So even if the air is stilled by a nearby boundary, the drag of the air on the free surface of the water will 
be negligible (provided the boundary is not too close). In the absence of a rigid boundary in the air, 
negligible error is made by treating the air as inviscid. 

Time In! 


Problems with Non -Zero Inertia 

L&L list only three flow problems in which both including viscous and inertial terms are important 
and in which exact solutions are known: 

1) rotating disk 

2) converging (or diverging) flow between nonparallel planes 
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3) the submerged jet 


Rotating Disk * 

We will now examine the solution of the first problem — 
the rotating disk — because it is used as a model system for 
transport experiments. An infinite plane disk immersed in a 
viscous fluid rotates uniformly about its axis. Determine the 
motion of the fluid caused by this motion of the disk. This 
problem was first solved by von Karmen (1921) using 
cylindrical coordinates with the /-axis coinciding with the axis 
of rotation. 


I 



I 


Define: 



where v = p/p [=] cm 2 /s 


is called the kinematic viscosity. 


Then: v r (r,z ) = rcoF(Q 

v e (r,z) = rcoG(Q 
v z (z) =a/v to H(Q 


p(z) = pwP(Q 

Continuity and N-S become: 

2F + H'=0 

r. F 2 +F'H-G 2 -F" = 0 

0: 2FG+HG'-G” = 0 

z: P'+HH'-H" = 0 



* See S:p93 (6th Ed); L&L:p79. 
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where the prime Q denotes differentiation with respect to Boundary conditions take the form: 


z=0: 

F=H=P= 0, G=\ 

z=°°: 

F=G=0 


An important property of this solution is that the z-component of velocity depends only on z. For this 
reason, convective heat and mass transfer problems becomes one-dimensional. 

This is perhaps the only problem in which there is flow normal to a wall where the mass-transfer 
coefficient can be determined analytically. For this reason, the rotating disk is a favorite tool of 
researchers in mass transfer. 


Creeping Flow Approximation 

Cone-and-Plate Viscometer 

The cone-and-plate viscometer consists of a flat plate, on which is 
placed a pool of liquid whose viscosity is to be measured, and an 
inverted cone which is lowered into the liquid until its apex just 
touches the plate. The cone is then rotated at some angular velocity, 
and the torque required to turn the cone or to keep the plate 
stationary is measured. 

Spherical coordinates are most convenient to describe this problem, 
since the surface of the cone and of the plate can be defined as 
0=const surfaces: 



surface of cone: 0 = a 

surface of plate: 0 = jt/2 

The cone is undergoing solid-body rotation (see HWK #4, Prob. 3): 
for 0<a : v(r) = flxr 

In spherical coordinates, the position vector is 
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r = re r (0,4>) 

while the angular velocity is (see figure and p49): 

Q. = £le z = £2[(cos0)e,. - (sin0)e e ] (56) 

Substituting into the expression for solid body rotation: 

0<a : v = Qxr = rCl sinO e<j, 

The principle direction of flow is the (^-direction. “No-slip” 
requires: 

at0=a: v^=rQ.sma 

v r = Vq = 0 

at 0=tc/2 : = v r = v 0 = 0 

The simplest velocity profile which is consistent with these boundary conditions is: 

^ = v^(r,0) 
v r = vq = 0 

Next, I’d like to argue that the pressure profile can be expected to be independent of <]). 

p =p(r,0) 

Since there exists only a (^-component °f velocity, fluid streamlines will turn out to be circles (the 
contour corresponding to r=const and 0=const). The circle corresponds to 0 < § < 2%. 

By analogy with the last problem, you might guess (incorrectly) that pressure must decrease along the 
direction of flow. However, in steady flow p cannot continously decrease with (|) for all (|). At the very 
least, pressure must be periodic in (f>; in other words, p(r.0.<])) = p(j\Q.<\>+2n). So any decreases in 
pressure over part of the cycle will have to be balanced by increases over the remaining part. Why 
should the pressure be higher at some points along the streamline than at other? There is no geometrical 
asymmetry with respect to cf> and no reason to expect any cf> -dependence in the pressure. 

The velocity profile automatically satisfies continuity: 

V • v = (rsin0) _1 3v ( |/9(|) = 0 
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Ignoring gravity, the Navier-Stokes equations in spherical coordinates become: 1 ’ 
r: -pv^/r = -dp/dr 

0 : -pvJ-cotQ/r = - r 1 dp/dQ 


i a f 2 ^6 


0 = px r — — + sin0 — — - 

r O "\ O "\ /~v ~\ /~\ 


sinO 50 


50 ) r 2 sin 2 0 


Notice that the pressure and velocity fields have been separated. We can first solve the ([(-component 
for the velocity profile and then substitute the result into the r and 0-components to solve for the 
pressure profile. Based on the boundary conditions, we might try a solution of the form: 

v <j)( r '0) = rf(Q) 

When this is substituted into the ^-equation above, the r-dependence cancels out, leaving a second 
order ordinary differential equation in/(0). The solution leads to: 


=r£2sina 


8 {a) 


where 


, . n 1, ( 1 + cos0 3 . 

g(0) = cotO +— In sin0 

2 ll-cosOJ 


This function is plotted in the figure at right. Notice that in the 
center (i.e. for 0 near n/2) the function is nearly linear. A 
more careful asymptotic analysis would reveal that 


lim{g(0)} = 2£ + o(e 5 ) 

it 1 ' ' 


where 


8 - 2 _e 


Notice that this is a linear function of 0 as expected from the plot above. 


* http://www.andrew.cmu.edu/course/06-703/NSE sph.pdf 
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Velocity Profile for Shallow Cones. Most cone- 
and-plate viscometers are designed with cone angles 
near to %I2. The reason for this will be apparent in a 
moment. Note the arc lengths re and re 1 on the figure 
at right. For shallow cones (i.e. as a—>%/2), re 
asymptotically becomes the vertical distance from some 
arbitrary point (r,0) in the fluid to the plate, 




while re \ becomes the vertical distance between the plate at the cone: 


lim {r£ 1 }=/i(r) 




So in this lim it we can use (61) to replace 


as a— >n/2: g(0) — » 2e = 2— 

r 

h 

and g(a) ->2e 1 =2— and sina — >1 

r 

Then (60) simplifies to linear shear flow (at least locally) 

z 

for a — >7c/2: va = rQ— 

' h 


cone 



plane 


Notice that the rate of strain is independent of position: 


for a— »Jt/2: 


rQ. rQ Q 

— — = -pr = = — = const. 

dz h{r) r£| £, 


This is an important advantage for a rheometer, since all the fluid experiences the same strain rate. 
Later, we will show that stress is also spatially uniform for shallow cones. 
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Torque. To use this device as a viscometer, we need to interpret torque 
measurements. So let’s try to evaluate the torque for a given velocity profile. 
Recall that torque is force times lever arm: 


T = (rsina)F 

Using vector notation: 

T = rxF 





top V LOW 


In our problem, the force F is not 
concentrated at one point, but instead is 
distributed over the surface of the plate. 

Let’s consider the contribution to force and torque from some 
differential element of surface having area da. 

Once the velocity profile is known, we can evaluate the stress field 
from Newton’s law of viscosity. Given the stress field, we can 
calculate the force on any differential solid surface element of area 
da from: 

dF = n • Tda 

where n is a unit normal pointing out of the body dF acts on. 
Similarly, we calculated the contribution to the torque by crossing 
this force with the local lever arm: 

d T = rxc/F = rx(n ■ I da) 


Let’s calculate the torque exerted by the fluid acting on the 
stationary plate. Then we want to choose n to point out of the plate or 


n = e. 


The net force exerted on the lower surface is: 


F = \ A dF = -| A e e • I da 


Similarly, the torque exerted on the lower surface is: 


T = | 4 rxc/F = j 4 rx(-e e • T)da 


(62) 


First let’s consider: 


e e ■ I - T&r e r + / B 0 e e + ^ 0(^0 
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To evaluate the torque, we need to first cross this vector with the lever arm, which is the position vector. 
In spherical coordinates, the position vector is: 


r (r,e,<|>) = re ; .(0,< )>) 

ix(e 0 - T) = r r 0r (e,xe r ) + rr 00 (e,xe 0 ) + /•T e(|) (e r xe (|) ) 


- r ^ee e (|) - r yr 00 e e 


Anticipating that the torque vector T will have the same direction as the angular velocity O, we dot both 
sides by e z which equals -e 0 on the plate (0 = n/2): 


z 

da 


= - e 0-[ rX (- e e-T)] = e e -[rx(e0-T)] = e0-[rr ee e^ -r7^e e ]= -/T e 


Now we substitute Newton’s law of viscosity:* 

7fl,h = Xfl, 



sin0 — 

( ,, V 

v 9 

+ 

1 

3v e 

r 

30 

sinO 

( ) 


sinO 

3c|) 


(63) 


What’s left is to substitute our velocity profile from (60) into (63), which yields 


C sin a 

Tqo = — — where C = -2pX2— — (64) 

sin n ° 1 n 1 

and where a is the cone angle. 

Q=n/2) will be independent of 
leaving: 

T z = hfi-rTQ^Qnrdr) 


On the surface of the flat plate, 0 = rc/2 and sin0 =1. T 0O (evaluated at 
4> ; then we can choose da to be a ring of radius r and thickness dr. 


T 




sina 

$(«) 


(65) 


where R is the radius of the region of the plate which is wetted by the fluid. 

Shear Stress and Torque for Shallow Cones. If the cone angle a is very shallow (i.e. a—>n/2) then 
T 0O will be practically independent of position: from (64): 


* http://www.andrew.cmu.edu/course/06-703/NLV sph.pdf 
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for any 0. This means the entire fluid experiences the same stress. This has a number of important 
advantages for rheological studies making the cone-and-plate viscometer one of the important 
viscometric flows. Substituting (61) into (65) 

lim [T l = f 7t|afl 3 — 

£->0 j £ 

Comment on Solution. There is a problem with our solution to NSE: when (60) is substituted 
back into (57) and (58), there is no single function p(r,Q ) which will satisfy them. In other words, (60) 
is not an exact solution . It turns out that (60) is a reasonably good approximation if Q. is not too large. 
The exact solution has the form: 

v = rflsina — 

S 

where 0(Zl 2 ) means that this term vanishes like <T> 2 as Q — >0. 

The centripetal force on fluid elements undergoing a circular 
orbit causes those fluid elements to be “thrown outward” in 
the +r-direction. Since fluid elements near the rotating cone 
are rotating faster than fluid elements near the stationary 
plate, we have outflow near the cone supplied by inflow near 
the plate. The resulting r- and 0-components of velocity are 
called secondary flow, whereas the original (^-component is 
call the primary flow. 

£l- (the secondary flow) vanishes faster than Q (the primary flow), so for small Q, the leading term 
is approximately correct. This is called the creeping-flow approximation. 



i_L 

(a) 


e,K +0 


( q2 ) 


Creeping Flow Around a Sphere (Re-»0) 

Let’s return to the problem of flow around a sphere 
(or motion of a sphere through a stagnant fluid). For jtf/k 
boundary conditions, we impose “no slip” on the surface 
of the sphere and far from the sphere the flow is — ► 
undisurbed: ^ 

at r=a: v = 0 

as r— >° o: v — > U 
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Here an exact solution to the Navier-Stokes equations is not possible. Of course, the vector equation 
can be converted into a scalar equations using the stream function, but that yields a 4th order nonlinear 
P.D.E. Although this could be solved numerically, considerable simplification can be obtained if either 
the viscous terms or the inertial terms can be neglected — even if they are not identically zero. One 
limiting case is creeping flow which corresponds to the lim it in which the Reynolds number is small (i.e. 
Re—>0). In this lim it the inertial terms in the Navier Stokes equations can usually be neglected. 


Scaling 

To show that inertial terms are neglibible, let’s try to estimate the order of magnitude of viscous and 
inertial terms for uniform flow at speed U over a sphere of radius R. 

1 V 2^ V = -Vp + ^|v 

inertia viscous 

We will use a technique called scaling (akin to dimensional analysis). We start by listing all the 
parameters in the problem. In this problem, the parameters are 

parameters: U. R, p and p 

A characteristic value for each term in the equations of motion is then written as a product of these 
parameters raised to some power: 


each term °c U a R b p c p cl 

For example, throughout most of the region, the fluid velocity is undisturbed: 

M~U 

where the symbol should be read as “scales like” In “scaling” we ignore any position dependence 
as well as any numerical coefficients, so |v| scales as U, although |v| might be significantly less than U 
near the surface of the sphere. From the boundary conditions, v changes from 0 at r=R to U at r=oo. 
To estimate the magnitude of the gradient Vv we need to estimate over what distance most of this 
change occurs. In the solution of the potential flow problem, the velocity profile along the rear 
stagnation line is 


v r (r, 0 = 0) = U 

We expect something similar for 
scalar components of Vv are 


( R\ 


COS0 


\r J J 

Stokes flow. Two of the nine 
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c)v f . 
dr 


D 3 

3 U cosO 

4 


u_ 

R 


1 dv r 
r 30 



3 

4 


sinO 


U_ 

R 


We again ignore the position dependence: we scale r as R and treat cos0 and sin0 as “1”. Then both 
components above scale as U/R as do other components of Vv associated with derivatives of v e . So 



Av _ U_ 
A r ~ R 


Lik ewise the velocity gradient can be expected to decay from a max im um value of U/R near the sphere 
surface to zero in the bulk; this change occurs over a distance on the order of R. 


|V 2 v 


A|Vv| _ % _ U 


A r 


R R- 


With these estimates, we can further estimate the magnitude of viscous and inertial forces: 


inertia = p|v- Vv| ~ p(f/) 




piu 

R 


Ivr2 I U 

VISCOUS = (l|V v| ~ |l 

R~ 

inertia _ p U 2 /r _ p UR _ ^ 
viscous iiu/R 2 

where Re is the Reynolds number. So as Re—>(). we should be able to neglect inertial forces 

for Re« 1 : 0 = -Vp + (iV 2 v (66) 

V • v = 0 

at r=R : v = 0 

as r— >°° : v U, p -> 

Eq. (66), called Stokes Equation, is common to all low-Reynolds number problems; it’s not specific 
to the sphere. A common trick to reduce the number of unknowns is to take the curl of both sides of 
the equation: 
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0 = -VxVp + pVxV^v 

But VxVp = 0 

Moreover, using identity F.l, 

V 2 v = V fe )-V x(Vxv) = -curl 2 v 

o 

for incompressible flow. Then (67) becomes: 

curl 3 v = 0 (68) 


Velocity Profile 

We might try to seek a solution having the form of potential flow: 

v = V4> 

since Vxv = VxV(|) = 0 


(69) 


(68) is automatically satisfied. But we know that the potential flow solution for the sphere does not 
satisfy the no slip boundary condition. On the other hand, the boundary conditions are axisymmetric: 

V0= 0, d/d§ = 0 

so we might seek a solution using the stream function: 

vM) 


v = V x 


r sin0 


'4» 


Computing the curl in spherical coordinates using the tables (http ://w w w .andrew .cmu.edu/course/06- 
703/Vops sph.pdf ) (also see HWK #5, prob. 1): 


curl 


f e^\|/ A 


r sinO 


= v = 


1 d\\r 


1 3 \| / 


J 


l 2 4'2% r f 

V r V 0 


ee 


Taking the curl a second time, using the same tables: 


curl" 


e ()¥ 
r sinO 


Vxv = K 


J 


rsinO 


(-£>) 


(70) 
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where E- is a partial differential operator given by: 

—,2 3 2 \i/ sinO d ( 1 9\|/3 

0 r 2 r 2 dOysinOdOy 

Since we have to take the curl twice more, it might look like we have a lot of algebra to look forward 
to. But, as it turns out, the rest is easy. We need to evaluate: 


curl 3 v ■ 


curE (V x v) = curl 


2 < 


K4 . 

rsinO 


where the second equahty is (70). The argument of this final curl-squared has exactly the same form as 
that of the left-hand side of Eq. (70), except that the scalar +\|/, which is a function of r,0, is replaced by 
-E 2 \\i. which is also a scalar function of r,0. So all I have to do is to replace \\i by -E 2 \\i on the right- 
hand side of (70): 


3 2 

curl v = curl 


e +(- £2 y) . 

rsinO 


®(j> 

rsinO 




rsinO 



To satisfy (68), which represents the curl of the Navier-Stokes equation, we choose the streamfunction 
to satisfy: 


E 2 (E 2 \y) = 0 

Translating the boundary conditions at r=R in terms of stream function: 

v r =0 : O\(//O0 = 0 at r=R 

v e =0: d\\r/dr= 0 

Translating the boundary conditions at r— >°o in terms of stream function (see HWK #5, Prob. la): 

as r— >oo ; \p — > ( l/2)t/r 2 sin 2 0 (71) 

The trivial solution V|t = 0 satisfies the P.D.E. and the b.c. at r=R. but not the b.c. at r— >° o. Based on 
the form of this nonhomogeneous b.c., we guess the solution has the following form: 

Try: x|/(r,0) =/(r)sin 2 0 

E 2 \\i = Zf2[/(r)sin 2 0J = ... = (/'”-2r 2 /)sin 2 0 = g(r)sin 2 0 

where: g(r ) =f"-2r 2 f 

Then: E 2 (E 2 \\f) = (g"-2r 2 g)sin 2 0 
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The Navier-Stokes equation become 


£ 2 (£ 2 y)= 0 : 

g"-2r 2 g = 0 

(72) 

def n of g : 

f"-2r 2 f=g(r) 

(73) 

b.c. ’s: 

f — > (1/2) Ur 2 as r — 0 

(74) 


./ =/’ = 0 at r=R 

(75) 


These two coupled second-order O.D.E.’s can be combined to produce a single fourth-order O.D.E. 
in f(r). The result is a Cauchy-Euler equation whose general solution is (see footnote on page 52): 

f(r) = c | r 1 + c 2 r + c 3 r 2 + C 41 4 

Applying b.c. (74): C 4 = 0 

c 3 = (1/2)6/ 

As r becomes large, terms which are proportional to higher power of r dominate those of lower power. 
To have the third term win over the fourth, requires us to kill the fourth term by setting its coefficient to 
zero. The remaining constants are evaluated in a straightforward manner. The result is: 


21 R 2 r 1 r 
N '(r,B) = UR- i — |- + i - 


4 R 2 \R 


r • 2 Q 
— sin 0 


v r (r,Q) = U cos0 

R f R 3^ 

v e (r, 0 ) = - 6 / l-j j — sin 0 

4 r 4 yrj 

The figure at right compares the streamlines for 
Stokes flow with those for potential flow. The 
streamlines correspond to values of the 


Stokes flow potential flow 



streamfunction which are uniformly spaced at 
about the same interval for both profiles. 

Displacement of Distant Streamlines 

From the figure above, you can see that streamlines in Stokes flow are displaced away from the 
sphere - even those streamlines which are quite distant - especially compared with potential flow, in 
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which the distant streamlines become straight. It turns out that all streamlines are displaced away from 
the sphere in Stokes flow. 

Consider the streamline in the figure at right which 
corresponds to \\i = \|/ 0 . Far upstream, the 
coordinates of the streamline correspond to r— >° ° 
and 0— >n such that rsinG — > const, which we will 
denote as y (see Hwk #5, Prob. la). The 
relationship between y and \|/ 0 can be deduced from 
(71): 



lim \|/(r,0) — » 

r— >oo ^ ^ 

e — >7i y 


¥ 0 =\ u y 2 


i R 

3 r 90 

, i 

( V \ 

r 90 

2“ 

l 

o 

4 R 

h 2 

\r J 



At the equatorial plane 0 = k/2, the r coordinate of the streamhne must satisfy (76): 

¥(^90>f)=¥o =UR 2 
Flim inating \|/ 0 between (77) and (78): 


2 2 3 1 R J 

y - r 90 ~T Rr 90 + T 

2 2 r 90 


or 


y = r 9o- 


i- 


(— 

v r 90 y 


+ ' 


1 


R V 

l r 90y 


(77) 


(78) 


For distant streamlines, r 90 will be very large compared to R , so that R/rg q « 1- Then the square-root 
in the expression above is unity plus a small correction, which can be estimated as the first term in a 

Taylor series expansion of the square-root function: V 1 + £ = 1 + ^-e+K 


y = t 90 


1- 


_R_' 

V r 90y 


+ 0 


_R_' 

V r 90y 


r 90 --i? 


where the “OfR/n^) 2 ” means that these terms decay to zero like (R/rgg) 2 as R/rgQ—>0. When Rh'gg is 
sufficiently small, we can neglect these terms compared to the others. This is how the second 
(approximate) equality above was obtained. 
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Finally lim {r 90 - y} = 6 = jR 

r 90 °° 

Thus, even streamlines infinitely far away are displaced a distance equal to % R. 


Pressure Profile 


Once the velocity profile is known, the pressure can be determined from 


Substituting v: 


Vp = pV^v 


dp 

dr 


3 \lURr 3 cos0 


-^t = lpURr~ 3 sin0 
r 30 2 r 


and then integrating with p—>p Xl at r— >°° as the b.c. 
yields: 


/ 3 r cosO 

p(r,Q) = p oa --[lUR—- U 

2 r ► 

This profile corresponds to a higher pressure on the moderate p 
upstream side of the sphere (0=71) than on the 
downstream side (0=0). Thus it appears as if a 
drag will be produced by this solution. 



-\ A pnda = ... = 2%\iRUk 


is called the form drag which arises from normal stresses. This however is not the total drag. More 
generally, we expect the net force exerted on the particle by the fluid is: 

F = j>n-T 'da= -jnpda + jn-Tda =6;t|a7?f/k 
a 142 43 3)42 43 

27tgtf[/k 4:t|VK[/k 
(form drag) (skin friction) 


where A is the surface of the sphere r=R 

n = e ; 

From symmetry, we expect that this force will have only a z-componcnt (the direction of bulk motion); 
so let’s concentrate on finding that component: 
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F z = j A k • (e r - T)da 

Substituting T = -pi + x 

e r - 1 = -p(e r ’l) + e r • x 
= -pe r + e r - x 

e,. • x = x, T e,, + x r0 e e + x^e^ 

For the determined velocity profile, Newton’s law of viscosity in spherical coordinates (W:146) yields 
XnjpO for all r and x, T =0 at r=R\ 

at r=R: e,.- T = -pe r + x r0 e e 

k • (e r • X) = -p(k • e,.) + x,. e (k • e 0 ) 

Using k = e, from (56): k • (e r • T) = -pcosO - x^sinO 

Finally, we can integrate using an azimuthal strip of width RdQ and of 
radius /?sin0: 

da = 2K(Rd\nQ)(RdQ) 

(79) becomes: 

n 

F z = 2kR~ |J-p(i?,9)sin0cos0- x r Q(/?,0)sin 2 oJr/0 
0 

Finally x r0 is expressed in terms of the velocity profile using Newton’s law of viscosity and subsequent 
integration yields: 

F z = 2%[i RU + 47ip RU 

We have already noted that the contribution from the pressure profile is call form drag. The second 
contribution is call skin friction . The total drag force is 

F drag = 67ip/?U (moving reference frame) 
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which is known as Stokes Law (1850).* In the expression above, U is the velocity the distant fluid 
flowing around a stationary sphere. If we switch back to the original (stationary) reference frame (see 
page 48f), Stokes Law becomes 


F drag = -67t[X/?U (stationary reference frame) 

where now U is the velocity of the sphere moving through otherwise stagnant 
fluid. Note that the drag force acts in a direction opposing the direction of 
motion of the sphere. 

Experimental results for the drag force around submerged bodies are 
usually expressed in terms of a dimensionless drag coefficient. The quantity that is used to make this 
drag force dimensionless is 




[=] 


K.E. 

vol 


force 

area 


This quantity also can be shown (according to Bernoulli’s equation) to represent the A p required to stop 
fluid which is flowing at speed U. Multiplying this by the projected area gives the force required to stop 
the flow, which would otherwise pass through the sphere: 

n Fdrag 

D 19 7 

jp U 2 x n 

projected area 


where nR 2 is the projected area of the sphere; in other words, kR- is the area of the sphere’s shadow 
cast along the direction of flow. Using this definition, Stokes equation for the drag force on a sphere 
yields: 


C D 


U_ 

Re 


where 


Re.H 


* Sir George G. Stokes (1819-1903), bom in Shreen, Ireland, educated at Cambridge; theoretical 
physicist. 
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Comparison with experimental results co nfir m that this works very well for RecO.l. The reproduction 
above (taken from BSL, pl92) uses the diameter for Reynolds number, rather than the radius. The 
“friction factor” on the y-axis differs by a factor of two from the drag coefficient defined above. 


Correcting for Inertial Terms 

For larger Re , Stokes law underestimates the drag force. Of course, this is due to the increasing 
importance of inertia which has been neglected. A number of investigators have attempted to extend the 
validity of Stokes law by including inertial terms in the analysis. We will now summarize some of these 
and hint at the difficulties involved. 

First, note that as Re — >0 for flow around a sphere, the NSE approximated by 


pV^v = Vp 


which is called Stokes equation. Stokes (1850) solved this equation to obtain: 

F = 671 pRU. 

The usual approximation to solving differential equations containing a small parameter (i.e. Re ) is to 
perform a perturbation expansion. Basically, the idea is to perturb the small parameter away from 
zero value by means of a Taylor series. 
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Perturbation Expansion 

time out : Let’s illustrate the main idea behind this powerful mathematical technique by means of a 
simple example (Example 25.1 from Greenberg). Find the solution y(x ;e ) to the following ODE 
involving a small parameter which we denote as e : 



y' + y + ey 2 = 0 

(80) 

subject to the initial condition: 

y(0) = cos e 

(81) 


where the independent variable spans the range 0 < x < oo and the parameter 8 is small: 0 < e « 1. 
Note that the solution of this problem is almost trivial in the special case of e=0 (first-order linear ODE) 

y(x ;0) = e~ x 

but obtaining a solution when e ^ 0 is more challenging (because the ODE becomes nonlinear). The 
general idea behind a perturbation expansion is to “perturb” the easily obtained solution away from e=0 
by seeking a solution in the form of a Taylor series expansion about e=0: 


y(^ e )= |(f>§) 

>’o W 


+ 


dy(.r;£) 

de 

1 41 2 4 e $P 

>’i ( x ) 


8 + 


1 d 2 y(v;e) 

2! f j £ 2 

1 4 4 2 4 #3° 

37 M 


8 2 +K 


Of course, this assumes that y(x,e) is “analytic” about 8=0. ^ It remains to be seen if the solution has 
this property or not. Renaming the unknown coefficients (i.e. the partial derivates), we look for a 
solution having the form 


y(x,e) = y 0 (x) + y^e + y 2 (x)e 2 + ... (82) 

Note that the coefficients y 0 , y\, y 2 ■■■ arc not functions of the parameter 8. Substituting (82) into (80): 

Cv o + + y 2 e2 + — ) + Cvo + y\ e + yi e2 + •••) + £ Cvo + .7i e + y 2 £2 + — ) 2 = o (83) 

Next we expand the squared sum by distributing each term of the first series over each term in the 
second series: 


Cvo + .7i e + y 2 £ 2 + •••) 2 - Cvo + Jt 8 + y 2 e2 + — )0 ; o + 1 7 ! 8 + yi^ 2 + •••) 


* “Analytic” means that these partial derivatives exist and that the series converges to y(,i\e). 
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Jo + Jt e+ J2 g2 +k 

Jo + ji £ + J2 £ ~~ +K 
Jo + JoJl £ + JoJ2 e ~ +K 
JoJl £ + Jf £ ~ + K 

JoJ2 e2 + K 

Jo +2joJl e +(j 2 + 2 Jo J2 j £ ~ +K 

Substituting this result for the squared sum in to (83) and collecting terms of like power in e : 
Solution continues . . . and will be completed in the next revision of the notes. 


First, let’s write the equation in dimensionless form: we will denote the dimensionless variables using an 
asterisk: 



V*= RV 


±- P~ P°° 

[lU/R 


£ = Re = 


p UR 
F 


The Navier-Stokes equations for steady flow become: 

£y*. ^ * y * = ^ y* — ^ * p * 

We then look for a solution having the form of a Taylor series expansion about 8=0: 

v*(r*,0;e) = v o (r*,0) + £ v^r*,©) + e 2 v 2 (r*,0) + ... 

Similarly for the pressure profile. Substituting this infinite series into the Navier-Stokes equation 
(dropping the *’s) 

£(v 0 + £v j +L y ( Vvq + eVvj +L ) = (V v 0 + £V ~ v ^ +L j — ( V pq + £Vpj +L ) 


0£° +(v 0 -Vv 0 )£ 1 +L =(v 2 v 0 -Vp 0 )£° +(v 2 v 1 -Vp 1 )£ 1 +L 


Next we bring all terms to the right-hand side of the equation. We then collect terms of like order in the 
small parameter (i.e. terms which are multiplied by £ raised to the same power). To obtain zero for the 
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sum for all values of e we cannot rely on cancellation of positive and negative terms.* Instead, the 
coefficient of each e n must vanish separately. This leaves the following equations: 

e°: 0 = V 2 v 0 - VpQ (84) 

e 1 : vp-Vvo = V 2 vj - V/^ (85) 

and so on ... Note that (84) is just Stokes equation. The solution for v 0 can then be substituted into 
(85) leaving a linear equation to be solved for V| and /q. 

This is the approach used by Whitehead (1889). v In many problems, this procedure works. 

Unfortunately, the solution for the higher order terms in the current problem cannot satisfy the boundary 
conditions. This result is known as Whitehead’s Paradox. ♦ Another method must be used. 

As an alternative, Oseen (1910) used an entirely different approach. He approximated the inertial 
terms and solved: 


pU- Vv = pV^y - Vp 


thus obtaining: F = 67ipMJ[l + (3/8)ite] 

In principle, one could refine the solution further by substituting the resulting solution for v in place of U 
and then re-solving for an improved v. In practice, although the Oseen equations are linear, their 
solution is sufficiently difficult that no second approximations are known. 

In 1957, Proudman & Pearson obtained the next order correction using a different technique 
called matched-asymptotic expansions . In this technique a different form for the expansion is sought 
near the sphere (which is called the “inner expansion”): 

r~R: x> = XQ l (r,Q)+Vi l (r,Q)(Re)+V 2 l ( r fi)(R e ) 2 +-- 


* If the sum of different functions of e vanished for one particular value of e (as a result of cancellation 
of positive and negative terms), then this same sum of functions evaluated at a different value of e would 
not necessarily vanish. The only way we can guarantee that the sum vanishes for every value of e is to 
make every term vanish for value of e . See also footnote on page 27 which concerns integrals rather 
than sums. 

v Al fred North Whitehead (1861-1947), English mathematician and philosopher, who collaborated with 
Bertrand Russell on Principia Mathematica (1910-13) and, from the mid-1920s, taught at Harvard 
University and developed a comprehensive metaphysical theory. 

♦ See Van Dyke, pi 52-3. 
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and far from the sphere (which is called the “outer expansion”): 

t»R: \° = Vo o (p,0)+v 1 o (p,0)(i?c)-t-V2 o (p,0)(i?c) 2 +... 

where p = rRe 

Hie inner expansion is identical with Whitehead’s which can be made to satisfy the no slip condition at 
r=R but the result does not have the correct form far from the sphere. Rather, the outer expansion is 
made to satisfy the boundary condition far from the sphere and then appropriately match with the inner 
solution to determine the remaining integration constants. The result is: 

F = 67ipRU[l + (3/8 )Re + (9/40)Re 2 lnRc + ...] 

The term inside square brackets can be thought of as a correction to Stokes equation: 

Re= 0.01 0.1 1.0 

[...]= 1.004 1.032 1.375 

Comparing the bracketed term with unity gives some idea of the error incurred by neglecting inertia. 


Flow Around Cylinder as Re-»0 


Now let’s look at the analogous problem of uniform flow normal to a cylinder at very low Reynolds 
number. If we drop the inertial terms in the Navier-Stokes equation, we obtain: 


pV^y = Vp 


V • v = 0 



r— >oo ; 


U 



r=R : v = 0 

For flow normal to the cylinder, we expect the velocity profile to 
correspond to 2-D flow: 



end view of 
cylinder 


v z = 0 and d/dz = 0 

So that a solution can be found using the stream function: 

v = Vx[\|/(r,0)e,] 

For slow enough flows, we expect that inertial terms can be neglected, leaving Stokes Equation (66). 
After the pressure is e li minated (by taking the curl of Stokes equation), we obtain a single equation in 
the unknown velocity profile : 
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recall (68): curl 3 v = 0 = curl 4 (\|/e z ) = V^(V^\\f)e z 

The boundary conditions are determined in a manner similar to those for a sphere (see Hwk #7, Prob. 


1): 


r —>°° : 

\|/ — » f/rsin0 

r=R: 

3v| t/dr = d\)//30 = 0 


Based on the b.c. at r—> °° (which is the only nonhomogeneous part of the problem), the stream function 
should have the form: 


\|/(r,0) =/(r)sin0 

Requiring V2(v2\j/) = 0 

generates an ODE for /(r) whose general solution can be obtained. Unfortunately, none of the 
particular solutions can satisfy all of the b.c.’s (see HWK #7, Prob. 1). This is known as: 

Stokes Paradox (1850) - Stokes equation for uniform flow normal to a cylinder has no 
solution. 

As it turns, the inertial terms dropped by Stokes are not entirely negligible - no matter how small Re is. 
Lamb (191 1)* obtained a solution for the circular cylinder as Re—>0 using Oseen’s approximation: 


pU- Vv = pV^y - V p 


V • v = 0 


Lamb’s solution for the drag force per unit length of cylinder is: 

F _ 47t \xU 

L i . Re 
j-- y- In — 

2 * 4 

where y = 0.577... is Euler’s constant. Notice that, un lik e Stokes’ solution for the sphere, Re appears 
explicitly in this result. No matter how small Re, this logarithm term in the denominator is never 
negligible. 


* Sir Horace Lamb (1839-1934), English mathematician who contributed to the field of mathematical 
physics. He wrote Hydrodynamics (1895) which was, for many years, the standard work on 
hydrodynamics. His many papers, principally on applied mathematics, detailed his researches on wave 
propagation, electrical induction, earthquake tremors, and the theory of tides and waves 
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Boundary-Layer Approximation 


Flow Around Cylinder as Re^ °o 


= -Vp + 

inertia viscous 


We have just seen that restricting attention to the limiting case of very small Reynolds number allows 
an analytical solution to the Navier-Stokes equation by neglecting or approximating the inertial terms. 
Similarly, we might expect that, in the opposite lim it of very large Reynolds number, we might be able to 
obtain an approximate solution by neglecting or approximating the viscous terms. 


Let's return to the problem of flow normal to a cylinder, but at very large Re. If we just drop the 
viscous terms from the Navier-Stokes equation, we get Euler's equation for an ideal fluid. Recall that 
a solution which satisfies the differential equation and the boundary condition far from the cylinder is 
potential flow. From HWK #4, Prob. 4a: 


v r 


(r,Q) = U 




COS0 


v e 0%9) = —U 




sin0 


At r=R : v r (R, 0) = 0 

and v e (R,0) = -2f/sin0 


which violates the no-slip boundary condition. An exact 
solution of the NSE's can be obtained numerically. 
Comparing the solution for Re»l to the potential flow 
solution (see figure at right), we see that the exact solution 
follows potential flow everywhere except in a narrow region 
near the surface of the cylinder where the exact solution 
turns downward in order to satisfy the no-slip condition. 



In the potential flow solution the velocity gradient goes lik e 


outside b.l.: 


3v e U 
dr ~ R 
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This is also the behavior of the exact solution outside the boundary layer. By contrast, near the surface 
(inside the boundary layer) as the Reynolds number increases, the velocity gradient gets steeper. A 
closer analysis (which we will perform in a few lectures) reveals 

inside b.l.: lim {- 7 ^- } = — -y/Re 

Re— >00 ^ or r= R) R 

The thickness of this region in which the two solutions differ decreases as the Re gets larger. This region 
is known as the: 

boundary layer, a very thin region near to a boundary in which the solution has a gradient 
which is orders of magnitude larger than its characteristic value outside the 
region. 


Mathematical Nature of Boundary Layers 

Boundary layers arise in solutions of differential equations in which the highest order derivative is 
multiplied by a small parameter. To illustrate the mathematical singularity which results, consider a 
simple example: 

Example: find asymptotic solution to the following problem as e — > 0: 

ey" +y' + y =0 

subject to: y( 0 ) = 0 

y(l) = l 

Problem is to find the asymptotic behavior of y(x) as £— > 0. This problem was presented Prandtl* (the 
father of boundary-layer analysis) to his class on fluid mechanics at Goettingen U. during the winter 
semester of 1931/2. 

Solution : For £ sufficient small, you might guess that the first term can be neglected, leaving 


y' + y = 0 


whose general solution is: 


y(x) = Aexp(-x) 


* Ludwig Prandtl (1875-1953), German physicist who is considered to be the father of aerodynamics. 
His discovery (1904) of the boundary layer, which adjoins the surface of a body moving in air or water, 
led to an understanding of skin friction drag and of the way in which streamlining reduces the drag of 
airplane wings and other moving bodies. 
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Now A can be chosen to satisfy either of the boundary conditions, but not both. To satisfy y(0) = 0, we 
must choose A=0 which does not satisfy y( 1 ) = 1 : 

y(0) = 0 — > A=0 — > y(x) = 0 for all x — > y(l) = 0^1 

On the other hand, if we choose A-e to satisfy y(l) = 1, then we cannot satisfy y(0) = 0: 


y(l) = 1 — > A=e — > y(x) = exp(l-x) for all x — > y(0) = e * 0 

The reason we can’t satisfy both boundary conditions with this approximation is that, by neglecting the 
first term, the order of the differential equation reduced from 2 to 1. With a first order O.D.E., we can 
only satisfy one boundary condition. Thus £=0 is singularly different from £ being arbitrarily small, 
but not identically zero. 


£— 0 — > O.D.E. is 1st order 
£ — > 0 — > O.D.E. is 2nd order 

Now the exact solution to this linear O.D.E. with constant coefficients is easily determined: 


y(x; £) = exp 


1 - x 
~2e 


sinh 

x r — — 
— V 1 — 4e 
Ue ) 


( 1 3 


2e 




Comparing the exact solution to that obtained by neglecting the term containing the small parameter for 
the case of £=0.05, we see that the approximation is good except near x=0. For smaller £, the region 
in which the exact and approximate solutions differ shrin ks . To see what's happening to cause this 
problem, let's take advantage of knowning the exact solution and deduce its asymptotic behavior inside 
the boundary layer. In particular, let's look at the initial slope: 
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/(0;e) 


v/l -4e 


exp 


m 

v 2 £. 


2£ 


sinh 


— Vl-4£ 

v 2e 


Now consider the limiting behavior as £— > 0. 


Vl-4£ -» 1 


and 


1/2e — » °° 


lim[V(0;£)] = lim 

e^0 e^0 


1 

2 £ 


( 1 3 


exp 


V 2£ j 


sinh 


m 

v 2 £. 


OO 

OO 


Qearly this lim it is highly indeterminant as both exp and sinh “blow up” fairly quickly as their arguments 
become large. This indeterminancy cannot be resolved with the help of L’Hopital’s rule because, no 
matter how many times you differentiate exp or sinh, they still “blow up.” The indeterminancy can be 
resolved instead by comparing the asymptotic behavior of sinh with that of exp: 


Recall the definition of sinh: 


so 


sinhz = 



exp(z) 2 
sinh(z) 


as z— > +°° 


as z — > +°° 


As z — > +°°, the second term of this definition becomes negligible compared to the first. Thus the ratio 
of exp to sinh approaches 2 and our expression above becomes: 

y'(0;e— > 0) « 1/e -> °° 

So the derivative of our function at this boundary is very strongly dependent on the value of the small 
parameter. In particular, the derivative is not bounded in the lim it e— > 0. This singularity is the essential 
nature of any boundary layer. In the above analysis, exp/sinh was bounded although both functions 
become unbounded as their arguments becomes large. We say that the singularity of these two 
functions in this limit is of the same order and we write this as: 

sinh(z) = 0(e~) as z— > °° 
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which is called the Big "Oh " notation (Bachman-Landau*). More generally, when we say that 

/(.*) = 0[g(.r)] as x— > a 

we mean f(x)/g(x) is bounded as .r— > a 

Likewise, we could say that: y'(0,e) = 0(8 _1 ) as e— > 0 

We could also show that y"(0,e) = 0(8 -2 ) ase->0 

which means that the second derivative blows up even faster the the first. This notation gives us a 
convenient way of describing the strength of a pole. 

at x=0: 8y" + + v = 0 (86) 

o(e _1 ) o(e _1 ) 0 

Thus at the boundary, the first term of the O.D.E. is the same order of magnitude as the second term. 
So no matter how small 8 becomes, the first term can never be neglected. This is the root of the 
problem. One way to obtain the correct order in e of the solution using a Taylor series expansion 
would be to transform the independent variable: 


Let 

X = x/e and Y(X\ e) = y(x; e) 

Then 

,, _dy _ dY clX _ i dY 
dx clX dx clX 


and 


// 


y 


d 2 y 
dx 2 


d | 

fdy 3 

d | 

( -t dY 3 
8 

d | 

( -1 dY ) 

dx 1 

v dx j 

dx ' 

l dX ) 

l ~~dx' 

v dX J 


clX 

dx 

-i 



d 2 Y 

clX 2 


Note that if dY/clX and d-Y/dX- are O(e 0 ) as e — ■> 0, we will obtain the correct order for y' and y". 
This is the basic idea behind “stretch transformation” which is part of the “inner expansion” which will 
presented in the second half of the next section. 


* This is the German mathematician Edmund Landau, not the Russian physist Lev Davidovich Landau. 
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Matched-Asymptotic Expansions* 

Matched-asymptotic expansions is a very general technique for coping with singularities lik e 
boundary layers. MAE is one type of singular perturbation expansion. Greenberg describes it as 
“one of the most important advances in applied mathematics in this century.” Although the beginnings 
go back to the 19th century and such names as Lindstedt and Poincare,* it was not until the 1960's that 
singular perturbation techniques became part of the standard analytical tools of engineers, scientists and 
mathematicians. Since boundary layers frequently arise in transport phenomena, let’s apply this 
technique to solve the simple problem posed by Prandtl. 

EXAMPLE: Use Matched Asymptotic Expansions to find the asymptotic behavior of the solution 
y(x;e) to the following problem as £— > 0. 


ey" +y' + y =0 


(87) 


subject to: y(0) = 0 

y(l) = 1 

Solution '. Following Prandtl (1905), we divide the domain into two regions: 

inner region: 0 < x < 8, y = y l satisfies inner b.c. 

outer region: 8 < x < 1 y =y° satisfies outer b.c. 

where 8 is the thickness of the boundary layer located near ,v=0. Within each region, we seek a 
solution which is a Taylor series expansion of the function y(x,e) about e=0: 


y(.v;e): 


?o(*) 


+ 


dy{x;e) 
de 

1 41 2 4 e $P 

yii x ) 


£ + 


1 d y(x;e) 

2! a £ 2 

1 4 4 2 4 # 3 ° 

37 M 


8 2 +K 


• The outer problem The solution outside the boundary layer can often be found as a regular 
perturbation expansion -. 


* The example problem introduced above was solved using MAE’s by Greenberg, p508ff. 

*(Jules) Henri Poincare (1854-1912), French mathematician, theoretical astronomer, and philosopher 
of science who influenced cosmogony, relativity, and topology and was a gifted interpreter of science to 
a wide public 
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y°(x,e) = y 0 (x) + V](x)£ + y 2 (x)e 2 + ... 
which is just a Taylor series expansion about 8=0. (88) into (87): 


e(yo" +e 3 ; l "+•••) + Cvo +e 3 ; i ,+ -") + Cvo +e 3 ; l + -") - 0 


Collecting terms of like power in 8 : 

O’o '+>!))£ 0 + (vo'+Vi '+3'i) e 1 + - = 0 

In order for this sum to vanish for all values of 8, each coefficient must separately vanish: 


0; 

y 0 ' + yo = ° 

(89) 

1. 

yf+yi = -Jo” 

(90) 


and similarly for higher order terms. Note that we have succeeded in obtaining a set of O.D.E.’s for the 
set of coefficient functions whose solution can be easily uncoupled. If we start with (89), VqCi') can be 
determined so that it is known when we solve (90). The outer solution must be subject to the outer 
boundary condition: 


y(l) = l 


Expanding this in a Taylor series about 8=0: 

y 0 (l) + yi (l)8 +v 2 (l)e 2 + ... = (l)e° + (0)8* + (0)e 2 + ... 

Thus y 0 (l) = l (91) 

yi(l) = 0 

and so on. The solution to (89) subject to (91) is: 

>’oW = exp(l-x) (92) 

We could now substitute (92) into (90) and obtain the solution for V|(v). Similarly, we could obtain 
y 2 (x), yfx) and so on. If we stop after the leading term, the outer solution becomes: 

y° = y 0 (x) + 0(8) 

y°(x,e) = exp(l-x) + 0(e) as e->() (93) 

• The inner problem To cope with the boundary-layer, we transform the independent variable 
using a stretch transformation , whose general form is 

X = 1 - X ° (n> 0) 
r.n 
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where xq is the location of the boundary at which the boundary layer arises. In this transformation the 
distance from the boundary Cy-a'q) is magnified (“stretched”) by an amount depending on the small 
parameter. In our problem, the boundary layer arises at a'q=0, so the transformation becomes 



Rewriting (87) in terms of the new variables X and Y(X) = y(x): 

/ _ dy__ d^dX_ _ 

dx dX dx 

Similarly y" = £~^ n $ 

where denotes dY/dX 

Substituting into (87): e 1_2,1 $4 Y = 0 (94) 

The purpose of the stretch transformation is to make all 0(8 °): 

as £ — >0 (X=const): Y, = 0(8°) 

whereas y,y',y" = O(8 0 ),O(8‘ n ),O(8 _2n ) 

Now we are in a position to choose a value for n. Recall from the exact solution [see (86)] that the 
second-derivative term is not negligible inside the boundary layer. This is generally the most important 
consideration in choosing n: we select the value of the parameter n such that we do not loose the term 
containing the highest order derivative (afterall, dropping the highest order derivative is what lead to the 
outer expansion, which we have shown fails in the inner region). To keep the highest order derivative, 

• This term must be lowest order* in 8 (otherwise it will be swamped by a lower order term) 

• This term must not be the only term which has this order (if this is the only term of that order, O.D.E. 
requires it to vanish identically) 

As a first attempt, we might try to make all terms of the same order: 

1 - 2/z = -n = 0 


*The “order” of a term with respect to a small parameter e (as opposed to the order of the derivative) 
refers to the exponent (power) to which e is raised. For example, we say that a term which tends to 
vanish like e n as 8 — >0 is “of order n.” 

Copyright © 2000 by Dennis C. Prieve 



06-703 


124 


Fall, 2000 


which is impossible for any single value of n. Next, we might try balancing first and third: 

1-2/7 = 0 or zz=l/2 

Orders (i.e. the exponent of 8) of the three terms then are: 

0 , - 1 / 2,0 

This is no good, because highest order derivative (i.e. the first term) is not lowest order in 8. Finally we 
try to balance first and second term: 

1-2/7 = -n, or zz=l 

Now the orders of each term is: - 1 , - 1 , 0 

which is OK. Using n= 1, (94) becomes (after multiplying by 8): 

eF = 0 (95) 

Instead of (88), we seek a solution in side the boundary layer which has the following form: 

for x<8 : y'(x,e) = Y(X,e) = Y 0 (X) + 7,(X)8 + Y 2 (X)e 2 + ... (96) 

(96) into (95) and collecting term of like order in 8 : 

8°: ^+1^ = 0 

whose general solution is: }q(A) = A + Bcxp(-X) 

Applying the inner boundary condition: Tq(0) = 0 

we can evaluate B- -A. This leaves us with 

y 0 (X) = A[l-exp(-A)] (97) 

Similarly, we could determine Y](X), Y 2 (X) and so on. If we stop after the leading term, we have for 
our inner solution: 

y' = A[l-exp(-X)] + 0(e) 

This remaining integration constant A must be chosen so as to match the inner and outer solutions. One 
possible choice is to take the outer lim it of the in ner solution [denoted (y‘)° ] and equate it with the inner 
lim it of the outer solution [denoted (y")'|: 
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lim |y*(X,£)l = lim |y°(.*;,£) 

1^2 4 43 J 1-44 2 4 



which is called the Primitive Matching Principle and was originally used by Prandtl. Taking the lim it 
of (93) as x— >0 and the limit of (97) as X— >°°: 


The inner expansion becomes: 
which means that for x<8 : 
whereas for x>8 : 


A = e 

y‘ = e[l-exp(-X)] + 0(£) 
y ~ e[l-exp(-x/£)] 
y ~ exp(l-x) 

v t? 


A convenient choice of 5 is where y l 
and y° intersect. Of course, this 

intersection point depends on £ : (>’ ! )° 

8 = 8(e) 

which is 0(£) in this problem. From 
the figure at right, neither y 1 nor y° is a 
good approximation to y in the vicinity 
of 8. However, in most transport 
problems, the quantity of greatest 
interest is dy/dx at x=0 and dy'klx 
does seem to match dy/dx at x=() quite 
well. If required, y i and y° can be 
blended together to obtain a single smooth function over the entire domain: 

y c = yi -|- yO _ (y*)o 

which is called the composite solution. For the present example, this is obtained by adding (92) and 
(97), expressing them in the same independent variable, and subtracting e\ 



c l-x , 
y = e +e\ 


{\— e A/e j-e = e{e 


—x/e ' 
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The plot at right compares y c (the dotted curve), and the 
exact y (the solid curve). Note that y c is a reasonably 
good approximation to y. The agreement gets much 
better as e — > 0. Better agreement could be obtained for 
any e by including the next order term in the expansions. 
Greenberg (p508f) obtains the second term in both the 
inner and outer expansions; the corresponding composite 
solution for e = 0.05 is virtually indistinguishable from the 
exact solution. 


MAE’s Applied to 2-D Flow Around 
Cylinder 



Let’s now try to apply MAE’s to solve the problem of flow 
around a cylinder at high Reynolds number. First, let’s write the 
equation in dimensionless form: we will denote the dimensionless 
variables using an asterisk: 


y* = 


y 

u 

*- P-Po 


r* = 21 

R 


V* = RV 


P U' 


£ = — = Re -1 
p UR 




end view «r 
e\ Linder 


This nondimensionalizing differs from our previous attempt which was for the opposite limit of small 
Reynolds number (see pi 12 of Notes). First, we have used pU 2 to nondimensionalize the pressure 
instead of \xU/R. This is because the disturbance to pressure caused by flow is proportional to pU 2 in 
the potential flow solution (see Hwk #1, Prob. 4). The second difference is that e is defined as the 
reciprocal of the Reynolds numbers, rather than the Reynolds number itself. This choice makes e a 
small parameter in the lim it Re— >°o. 


The Navier-Stokes equations for steady flow become: 



v *.v * v* = eV * 2 v*— V*p* 

(99) 

and 

< 

* 

< 

* 

ii 

o 

(100) 

Boundary conditions are 



as r*— >°°: 

v* — > e x and p* — > 0 

(101) 

at r*=l: 

v* = 0 

(102) 
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Expecting a boundary layer to arise at r*=l as Re— >°o (or as e— >0), we will now use the technique of 
MAE to solve it: 


Outer Expansion 

In the outer region, the regular perturbation expansion in powers of e uses r * as the position 
variable: 


v*(r*,0;e) = v o (r*,0) + e v^r*,©) + e 2 v 2 (r*,0) + ... 

p*(r*,0;e) = p q ( r*,0) + e p i(r*,0) + e 2 p 2 (r*,0) + ... 

Substituting this perturbation expansion into (99) through (101), collecting terms of like order in e, and 
setting their coefficients to zero, produces a series of well-posed mathematical problems for the 
coefficient functions. The first problem (the only one we will worry about in this analysis) is (dropping 
the *’s): 

£°: v 0 - Vv 0 = -Vp 0 (103) 

V • v 0 = 0 (104) 

The outer expansion is required to satisfy the outer boundary condition: 

as r— >° «: v 0 — > e T and p 0 — > 0 (105) 

Although we should not generally require the outer expansion to satisfy the inner b.c.’s, when we later 
match inner and outer expansions (see footnote on page 130), the outer expansion will have to satisfy: 

v,^ = 0 at r=l (106) 

Note that the viscous term does not appear in this result because the lowest order viscous term is 
Ole 1 ), whereas other terms are 0(8°). It turns out that (103) through (106) is the same problem we 
previous solved by potential flow: 


v 0 = V4> 


where is chosen to satisfy (104): V 2 9 = 0 

For potential flow, Vxv 0 =0 and (103) becomes: 

v ( po + T ’») =0 

which is just Bernoulli's equation for the pressure profile. After imposing the outer boundary conditions 
in (105) (and v ,.=0 at r= 1) we get the potential flow solution (see Hwk Set #4, Prob. 4). Eventually, 
we will need the inner lim it of the outer solution for matching: 
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v°(r = l,6) = 0 
ve(r = 1,0) = -2sin0 


p°(r = 1,0) 


1 -4 sin 2 0 
2 


(107) 

(108) 


Inner Expansion 

For the inner expansion, we use a “stretched” radial coordinate: 

F = or r* = 1 + 8 n Y (109) 


and rename the tangential coordinate X = 0 

U nlik e our simple example problem in the last section, the stretching parameter n in this problem is not 
an integer. More generally, the inner expansion should be a power series in e” rather that a power 
series in e itself: 


v5(r*,0;E) = M O (X,y)+£' ! M 1 (X,F) + £ 2n M 2 (X,F)+K 
v*(f*,0;£) = v 0 (X,F) + £ m v 1 (X ,F) + £ 2,1 v 2 (X,F)+K 
p*(r*,e;£) = p 0 (X,Y) + £ n Pl (X,Y) + £ 2n p 2 (X,Y )+ K 
In cyhndrical coordinates for 2-D flow, the equation of continuity (100) becomes: 


^ * . y* : 


1 3 ( r * v *) i 


r* dr* 


dva 

H — — = 0 

* 30 


We constmct the first term using (109) through (1 1 1): 

= ^1 + £ f)(v 0 + £ Vj +K ^ = vq + (Fvq + vj )e +K 


d(r*v*) 


r *v r 


dY d 


dVn -r, 


dr* dr* dY 


v o +(^ v o + v 1 )e"+K J--^-e 


+ 


3v 0 3v t 


vo+Y-*r + 


dY dY 


£°+K 


(HO) 

(HI) 

( 112 ) 
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In the analysis above, we needed to expand \l + e n Y) as a power series in e". This was 
accomplished using the Binomial Series, which will be quite useful in later problems as well: 


,, ,a 1 a(a-l) 2 a(a - l)(a - 2) 3 

(1 + x) = 1 + cvc + x + x +K 


2 ! 


3! 


which is just a Taylor series expansion about v=0. The Binomial Series is known to converge provided 

W<i. 


Similarly 


^’9 _ du 0 +£ « }kl \ [K 


30 dX 


dx 


Hie continuity equation becomes: 


dv n — /! 
— —8 + 

dY 


3vi du n 3 
0 dY dX 


8°+K = 0 


To satisfy “no slip” at the inner boundary, we must require: 

at F=0: u 0 = Mi =K = 0 and v 0 = Vi =K = 0 

Setting the coefficients of each term separately to zero: 

d\’r 


£-»: 


dY 


— = 0 for all Y 


v 0 (T) = const = 0 


(113) 
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which means that vq must be constant inside the boundary layer. To satisfy the no-slip condition (1 13), 
the constant must be zero. Then v 0 must vanish everywhere.* The next term of the continuity equation 
is more useful: 


£°: 


or 


\o 

o 


, dvi , du 0 
dY dX 


= 0 


duQ_ + dv 1 
dX dY 


= 0 


(114) 


Next, we will look at the principle component of the NSE, which will be the component tangent to the 
surface. For steady 2-D flow, the 0-component of (99) is (BSL, p85): 


NSE 0 : 


dv, 


dr 


e + At 


dv e v r v e 1 dp 

— -H-+_L- »- = — + e| 

r dO r r dO 


d 

'l d(n'e)' 

dr 

r dr 


1 d~va 2 dv, 
+ - - +• 


r 2 d 2 0 


30 


Next, we transform each term using (109) through (1 12) with vq = 0: 


NSE 0 : 


dllQ dllQ _ 

i ’M iV? 

4 °) 4 °) 


^o+K+e 1 - 2 "^-^ 

SX 14 2 


Ol 


(«°) 


o 




(115) 


As a general mle, we don't want to lose the highest order derivative inside the boundary layer, so we 
want this term to be lowest order in 8, but not the only term with this order. The largest inertial terms 
are 0(8°), so we choose n such that 1-2/7 = 0: 

1-2/7 = 0: /7=| (116) 


After collecting like-power terms, the r-component of (99) is: 


NSE r 



+ K =-r-H Sp 0 3p l 

142^ ^ 

o(„-'r-) o(c°) 


+ K 


*This result could be used to evaluate the undetermined integration constant in the outer expansion. 
Since v 0 represents the leading term in the inner expansion for v r , this result means that the outer lim it of 
the inner expansion for v r is zero, and this must match with the inner lim it of the outer expansion. 
Indeed, our expressions for the outer solution already make use of this result. 
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The lowest order term in this equation is Bpq/BY. Since there is no other term with this order, we must 
require that: 


£- 1/2 of NSE r : 


M) 

BY 


= 0 


which integrates to yield: 


Po = c ( x ) 


where the integration “constant” c(X) can be evaluated by matching the outer lim it (Y—> °o) of this inner 
solution with the inner lim it (r—> 1) of the outer solution (108): 


P 0 


c(X) 


l-4sin X 


(117) 


We still have two unknowns: w 0 and vp We can formulate two equations from the above: 

> 2 . 


£°of (115): 


and (114): 


u 0 XT 1 + v l — -^9- = 4sin X cos X 


BX 


BY dY 2 

Buq_ + Bv l 

BX BY 


0 


(118) 

(119) 


No slip requires: 

atF=0: u 0 = v t = 0 (120) 

Matching the outer lim it of the inner solution with the inner lim it of the outer solution requires: 

as X- w 0 — > - 2sinX (121) 


Boundary Layer Thickness 

Recall the definition for “boundary layer” from page 117: 

boundary layer : a very thin region near to a boundary in which the solution has a gradient which is 
orders of magnitude larger than its characteristic value outside the region. 
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So how thick is the boundary layer for 2-D flows like the 
one we just analyzed? A typical profile for the tangential 
component of velocity inside the boundary layer is sketched 
at right. On the scale of this sketch, the velocity appears to 

approach [v'q \ which is the inner lim it of the outer solution. 

Actually, the nearly horizontal dotted line isn’t quite 
horizontal: it has a small negative slope, but its slope is so 
small compared to the slope inside the boundary layer, that 
the outer solution appears to be flat on the scale of this 
drawing. 



We might define the thickness of the boundary layer as the 

distance we have to go away from the surface to reach the apparent plateau. From the geometry of this 
sketch this distance, which we denote as 8, is approximated using 



8 





( 122 ) 


where all of the quantities in this equation have units. Let’s try to “scale” (122) and solve for 8. Recall 
that the outer solution corresponds to potential flow. From HWK Set #4, Prob. 4: 


vg(r,e) = U 



sinG 


so that 



lim vq (r,0) = 2 U sin0 



In dimensionless quantities, the inner solution is given by (1 10). Substituting n = 1/2, we have 

vj(r*,0) = u 0 (X,Y) +e 1 / 2 u 1 (X,Y)+K 


dy 


Udv* Q 
R dr * 


U dY d 
Rdr^dY 
P ~Yi 


w 0 (X,F) + £^w 1 (X,F)+ K 

R 

O+k] 

dY 



Ia°) J 


Dropping any multiplicative constants, we obtain 
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H^ £ -y 2 u 

dy R 


Substituting this result into (122) and solving for the boundary-layer thickness 8: 


U_ 

8 


£ 


-%v. 


R 


or 


8~JeR 


R 

VRe 


Thus 8 is proportional to R and inversely proportional to the square-root of Reynolds number. The 
boundary layer gets ever thinner as the Reynolds numbers increases. This is hue for all boundary layers 
in 2-D flows. 


Prandtl’s B.L. Equations for 2-D Flows 

Let’s now summarize the mathematical problem which must be solved to obtain the velocity profile 
inside the boundary layer. The mathematical problem is represented by equations (118)-(121). For 
clarity, let’s rewrite these equations using variables having dimensions [recall (98), (1 10)-(1 12) and 
(116)]. For example, uq is given by (1 10) and (98) as 


* v e 
u o = v e = ~jj 


while Fis given by (109) and (116) as 


Y = 


r*~ 1 l '- R e -l /2 = 


d/2 


R 


p UR 

d 


d/2 


r-R 

R 


so 


dY = 


r pUR ^ 2 

d 


dr 

R 


d z u Q p R 2 d 2 r>e 

uR~u~dy 


Thus 


Making similar transformations to dimensional vaiiables of each term, then multiplying both sides of the 
equation by pL/“//? , (118) becomes 


"V0_dv0 
R 30 


+ Vy 


3vq 

dr 


d 2 

d— 

dr 


pu- . . . 

1 4sin OcosO 

R 


( 123 ) 
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(119) becomes 


1 dv Q | 3v r _ Q 
R 30 3r 


(124) 


(120) becomes 


v r = v e = 0 at r = R 


(121) becomes 


Vfi = -26/sin0 at r— > oo 


Notice that (123) and (124) are approximations to NSE e and continuity in cylindrical coordinates. 


Although the above equations for the velocity 
profile inside the boundary layer were derived for the 
specific geometry of a circular cylinder, it turns out that 
very similar equations are obtained for any 2-D flow, 
provided we express them in terms of a local 
Cartesian coordinate system (x,y). 

• x = arc length measured along the surface in the 
direction of flow 


cylinder 


• y = distance from the surface measured along a 
normal to the surface 

Hie more general equations for any 2-D flow are given by 


3 v x dv x 


x dx d y J " 


3 2 Ft dp 0 

dy 2 dx 


3v^ + 3v z 

dx dy 


(125) 


where I s the inner lim it of the pressure profile in potential flow, which is given by (1 17): 


,(*)= lim[// f (r,0)] 


For boundary conditions, we impose “no slip” at the wall: 


at y=0: 


v x (.r,0) = v v (jc,0) = 0 


Remaining integration constants are evaluated by matching the outer lim it of the inner solution with the 
inner lim it of the outer solution: 


as y— 


v xi x ’y)-^ u o{ x ) 


wher U 0 (x) is the inner lim it of the outer solution for the tangential component of velocity: 
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U 0 {x)= lim [v 0 PF (r,e)l 

r— > Rl J 


is the inner lim it of the potential flow solution. 

The the main difference from one geometry to another is a different potential flow solution applies to 
different geometries. From Bernoulli’s equation, we have 

Po ix) Unix) 

— + — = const 

P 2 

Differentiating with respect to x and rearranging: 


1 dpo _ j dU{ 
p dx 2 dx 




Substituting into (125): 


v *ar +v ’ 


dv r d"v. 

— v 

dy dy 2 


v x _ TJ au o 
"TT" — On — : — 


0 L ^ I{) = known/ (x) 


, 9v y 

dx dy 


(126) 


where v = p/p. These PDE’s are called Prandtl’s boundary layer equations. Appropriate b.c.’s 
include 


at y=0: 


v r = v y = 0 


as y— >°o: 


v_ r -> Uq(x) 
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Alternate Method: Prandtl’s 
Scaling Theory* 


Now we will repeat the analysis of boundary 
layers using Prandtl’s analysis, which is more 
intuitive. Consider uniform flow normal to a long 
cylinder. Of course this is 2-D flow: 


= 0, d/dz = 0 



At high Reynolds number, we expect the viscous 
terms to become negligible everywhere except 
inside the boundary layer. Dropping the viscous 
terms from the problem yields potential flow: 


r-R>8 : 


v = V(|) 


where 8 is the thickness of the boundary layer. 
Prandtl’s analysis of this problem assumed that: 



• outside the b.l., viscous « inertia, such that the velocity and pressure profiles are those obtained 
from potential flow. Recall the potential-flow solution from HWK 4, Prob. 3: 


v 


PF 

r 


u 


l- 


\r ) 


COS0 


for r-R > 8 



= -U 



V r J 


sin0 


• 8«R so that fluid elements inside b.l. don't “see” the curvature of the cylinder (we call this 
“ Postulate #i”). 

• inside the b.l., inertia ~ viscous (we call this “ Postulate #2”). 

Defining a local Cartesian reference frame, the continuity and Navier-Stokes equations become: 


continuity: 




= 0 


*Schlichting, 6th ed., pi 17-21. 
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x : 


y- 


BVy 

VX lh +Vy ' 


dy 


1 dp 

‘p^ + V 


f d 2 v x + 3\ 


dx~ dy ‘ 


d\’ 


dv. 


y y 

V ~ ^ + V v - 

x -\ y 


dx 


dy 


1 dp 

— + V 

pdy 


( ~\2 -\2 

d v d v 

“ + 

dx 2 


dy 


where v = p/p. Next, Prandtl estimated the order of magnitude of each term, hoping that some can be 
dropped. This is called scaling. In this estimate, we are not concerned about factors of 2 or 3, but 
with orders of magnitude: we will try to guess the asymptotic behavior of each term in terms of the 
characteristic physical parameters. In particular, in scaling we try to express each term in the form of 
the product of the physical parameters raised to some power: In this problem the physical parameters 
are R, U. v, and 8. Thus we will try to express each term as 

pa jjb v c g d 


Let’s start with the primary (x) component of velocity. Across the boundary layer, v x must vary 
between zero at the surface (no-slip) and 

at y=0 : v x = 0 


at y~8 : 


rPF = 2C/sin0 


which is the potential flow solution. The assumption here 
is that the outer edge of the b.l. corresponds to the inner 
edge of the potential flow solution. 


So 


v x aU 


dv r Av , 


dy Ay 8 



aV , 

dy 2 


v dy 
Ay 


0-8 8 2 


Now x is measured along the surface of the cylinder. 
Thinking of v as the arc length measured from the 
forward stagnation point, then: 

dx - -R dQ 


Integrating with v=0 at 0=7t : 
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x =R(n-Q) 


dv e PF dO TT J 1 
dO dx \ R 


d2 v x _ 

d 2 vg F f 

dx 2 

do 2 ( 


Similarly, — f = V" 

dx 2 dO 2 

Next, we look at the equation of continuity: 


a 2 = (2 t /si„e)f-lf .v 

dx) \ R) r 2 


^y _ <^x .. ^ 

dy O.y /? 

_ . , , , f 3v 2U cosO Ud 

Integrating across the b.l.: v v = — —dy ~ y ~ — 

' J dy 14?^B ’ R 

dv x /dx 

_ d ( 2f/ycos0 3 dO _ Ud 
dx 30 v R J dx R 2 

- 7 * 

Estimating the inertial and viscous terms in the v-component of N-S: 
inertia: v x dv x /dx ~ ( U)(U/R ) = U 2 /R 

Vydvjdy » (Ud/R)(U/d) = U 2 /R 
viscous: vd 2 v JOx 2 ~vU/R 2 

vd 2 v x /dy 2 ~ vU/d 2 

Now let's summarize our results as to the magnitude of each term in the x -component of the NSE: 

dv r d\’ y 1 dp d Vy d Vy 

x • v — — + v — — = — + v — + v — 

12% 12% 1/3' 12% i|% 

U 2 u 2 ? vU VU 


vd 2 v x /dx 2 » vU/R 2 


Since d«R (according to “Postulate #1”), the first viscous term must be negligible compared to the 
second: 


d«R: 


vU/R 2 « vU/6 2 
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d 2 v x /c)x~ « d 2 v x /dy 2 


So we can neglect this term. 

Ignoring the negligible term, we have established the following orders of magnitude for the terms in the 
A' -component of the Navier-Stokes equation: 


x : 



1 fy , y d 2 Vx 

i $ i J v 3 

? vU 

8 2 


(127) 


Now let's apply Prandtl’s second postulate: Inside the boundary layer, viscous and inertial terms are of 
the same magnitude: 

inertia ~ viscous 

U 2 /R ~ \>U/'S 2 


This allows us to estimate the thickness of the boundary layer: 

S 2 = vR/U = R 2 /(RU/v ) 

or 8 ~—IL= (128) 

VRe 

where Re = RUN 

is the Reynolds number. Note that this correctly predicts that the boundary layer gets thinner as 

8— » 0 as Re—> oo 

The remaining term in (127) which has not yet been estimated is the pressure gradient. We can obtain 
some idea of its magnitude by looking at the potential flow solution, in which the pressure profile is given 
by Bernoulli's equation: 

for y>8 : p/p + v 2 /2 = const 

Now the kin etic energy can be decomposed into contributions from the x- and y-componcnts, whose 
orders of magnitude we have already estimated: 

V 2 = Vjc 2 + Vy 2 - V, 2 

Now v 2 ~U 2 while v v 2 ~( Ub/R) 2 , so v v 2 «v Y 2 which leaves: 
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Differentiating: 


p/p ~ const - v x 2 /2 
(1/p )dpldx ~ vJdvJdx ~ U 2 /R 


(129) 


So this means we cannot neglect the pressure gradient in (127) since it is the same order of magnitude 
as the inertial terms. Scaling of the terms of the v-componcnt of the Navier-Stokes equation yields: 


y- 


dV y 

12 % 

U- 8 
R 2 


■ + V 

1 


dVy 

v 2¥ 

Ul 5 
R 2 


y 


1 dp -|- v ^ 

l %3 

? vt/S 

R 3 


+ V 




, 2 

1 2 % 

vU 

8R 


(130) 


Substituting Eq. (128) for 8 shows that viscous and inertial terms again have the same magnitude: 


inertia • 


viscous 


>2 V U 


R- 


R- 


so 


vU vU U 3/ _ 3/ I / 

~ ~ u ' 2 R /-\)' 2 

8R R \\R 


media ~ viscous 


What remains to be determined is the pressure gradient. At most, the pressure gradient in this equation 
has the same order of magnitude as the other terms: 


(1/p )dp/dy < U-6/R 2 


(131) 


Comparing this with the pariial derivative in the other direction, by dividing (13 1) by (129): 

( dp/dy)/(dp/dx ) < '6/R — > 0 as Re—> oo 

which means that variations in pressure across the boundary layer are becoming insignificant at large Re. 
In other words, a good approximation would be to take: 

p=p(x) 

Since the pressure at the outer edge of the b.l. must correspond to that just outside, where potential 
flow occurs, we can calculate this pressure using the Bernoulli's equation and the potential flow solution: 

for y>8: p/\ p + v 2 /2 = const (132) 

In potential flow at r=R : y = 0 

v 2 (R,&) = y e 2 (/?,0) = U 2 (x) (133) 
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(133) into (132): 

p/p + U 0 2 / 2 = const 



( 1 /p )Bp/Bx + U 0 dU 0 /dx = 0 

(134) 

(134) into (127): 

dv r dv r 3~v r TT dUn 

v r — — + v v — — - v — = U a — = known 

x dx y dy Bx 2 dx 

(135) 


The right-hand side of this equation is known because Uq is the inner lim it of the potential-flow solution: 

U ° ( x )-^[ v e PF ( r ’ 0 )] 

Together with the continuity equation, we now have two equations in two unknowns: 

2 unknowns: v x and v y 

2 equations: (135), Continuity 

which are known as Prandtl's Boundary -Layer Equations for 2-D flows. 

Solution for a Flat Plate 

Reference : Schlichting, 6th ed., pl25-33, Whitaker 
p430-440. 

Before we continue with the analysis of flow 
around a circular cylinder, let's look at the simpler 
problem of flow tangent to a semi-i nf inite flat plate. 

The analysis begins by computing the potential-flow 
solution. 

• Step 1 : find potential flow solution 
If the plate is infinitesimally thin, the uniform velocity profile is not disturbed: 

P.F. : v = U for all (.v,y) 

and P =Po o for all (x,y) 

Of course, this doesn't satisfy “no slip,” on the plate, but then neither did potential flow around a 
cylinder. 

Step 2 : apply Prandtl’s b.l. equations 

From this potential-flow solution, we can calculate the U 0 appearing in Prandtl's equation: 
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U Q (x ) = v x PF (x,0 ) = U (a const.) 

U 0 dU 0 /dx = 0 

Eq. (135) becomes: v x dv x /dx + v y dvjdy = vd 2 v x /dy 2 

dvjdx + dvy/dy = 0 

Appropriate boundary conditions are: 

no slip: v x = v y = 0 at y=0, x>0 

and outside the boundary layer, we obtain potential flow: 

v x — > U as y — > oo 

Step 3 : re-write the b.l. equations in terms of the streamfunction 
We can contract the two equations into one by using the stream function: 

v = Vx[y(x,y)k] 

which automatically satisfies continuity. Written in terms of the stream function, the problem becomes: 

V/l fxy - Wyy = ^Vyyy (136) 

at y=0,x>() : \[q. = \\r y = 0 

as oo : \|/ -> Uy 

where the subscripts denote partial differentiation. 
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The development of the laminar boundary layer along a flat plate is visualized by the hydrogen 
bubble method. A fine electrode wire is introduced upstream of the flat plate and a voltage pulse is 
applied reapeatedly at regular intervals. The boundary layer thickness is seen to increase with 
the distance downstream for the leading edge. Photo taken from Visualized Flow , Pergamon, 

New York, pl7 (1988). 

Flow visualization studies show that the boundary- 
layer grows in thickness as you move downstream 
from the leading edge (see sketch at right). This 
can be rationalized by considering the inverse 
problem of a plate moving through a stagnant fluid. 

Focus your attention on an intially stationary fluid 
element in the path of the moving plate. Hie longer 
the fluid is disturbed by the moving plate, the more time there is for momentum to diffuse away from the 
plate. 



Time Out: Flow Next to Suddenly Accelerated Plate 


First consider the much simpler problem of a infinite plate 
suddenly put in motion at 1=0. Initially, the fluid and wall are 
at rest; but at time 1=0 the wall is set in motion in the x- 
direction with a steady speed U. The in itial and boundary 
conditions are: 

at 1=0 for all y>0: v x = 0 



at y=0 for all />(): 


v x =U 


as y— > 


v v — > 0 
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Hie solution can be expected to have the form: v x = v x (y,t), v y = v z = 0. The NSE becomes: 

dv v d 2 v r 


-2L = v - 


dt 


dr 


where v = p/p. Hie solution to this problem is well known: 5 ' 


v x {y,t) = Ue rfc 


' y A 

V4v7 


h 


where 


erfcT) = 1 — -/= [ 

. I <TT • 


.... r s11 ‘ 

Jn A i [o as r| 

1 4 ( 4 2 4 43 

erf r) 


is the complementary error function', the integral itself is the error function. 

At any fixed t, the velocity decays monotonically from U at y=0 to zero as y— >°°. As t gets larger, 
more and more fluid begins to move; we say “the motion penetrates deeper into the fluid.” Suppose we 
wanted to know how far from the wall we have to go before the velocity drops to 1% of the wall value. 


: (y = 5,f) = O.Olf/ : 


(M 


( 


u 


= erfc 


A 


erf 


a/4v7 I 


V 


>/4v? 

1 - erfc 


= 0.01 


\ v -r V ' J 

( 

V 


-4= 1 = 0.99 
V4 w J 


Looking up the appropriate value in a table of error functions, we obtain: 

= 1.825 


or 


v/4vt 
8 = 3.65 Vw 


(137) 


This gives us some idea how far the motion “penetrates” into the stagnant fluid. Notice that the 
penetration depth increases with the square-root of time. 


see BSL, pl24f. 
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Time In: Boundary Layer on Flat Plate 


Now let’s apply this result to our original problem: the boundary-layer next to a semi-i nf inite flat 
plate. Fluid farther downstream from the leading edge has been in contact with the moving plate longer; 
and estimate of the “contact time” is: 


t = x/U 


(138) into (137): 


s °c ^^vt ,/ v '/ 


'U 


(138) 

(139) 


• Step 4 : solve b.l. equation using “similarity 

transform” 


U 


a 

a 

=a 

=a 



31 





6{x) 


Judging from the boundary conditions alone, we T 
would expect v x to vary from 0 at y=0 to U at y=8 — ► 

with sort of a parabolic shape, with 8 getting larger ► 

as we move downstream. Notice that the basic 
shape of this profile is not really changing with x, 

only the range of y- values over which the solution departs from potential flow is increasing as we move 
downstream. This suggests a solution of the form: 


v x IU =/' (y/8) 

Let's define a new independent variable which is scaled to the boundary-layer thickness: 


(140) 


Let 


S(.r) -Jvx/U 


= r|(.r,y) 


To get this guess in terms of the stream function, recall: 

5\|/ 

dy 


v * = 


x,y 


thus 


h 

¥ = J % ^z§') f = Ub \ f 'I 1 !) dr \ = U8 f I 1 !) 

x,y=0 jjf 8 dr) 0 


substituting 8 from (139): \|/(x, y) = -JvUxf(x\) 

In terms of/(ri), the boundary-layer equations (136) become: 

ff " + 2/ '"= 0 


( 141 ) 


subject to: 


/=/' = 0atri=0 
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/ ' — > 1 as r| — > oo 

Now we have an ordinary two-point boundary-value problem, which can be easily solved. This 
solution was first obtained by Blasius (1908), who was a Ph.D. student of Prandtl. 



Notice that the velocity component normal to the plate (i.e. v y ) does not vanish far from the plate: 



v y (*.y) 

U 


= 0.8604. 



* 


0 


Consider a fluid balance around the shaded 
rectangle in the figure at right. Owing to the need 
to meet the “no slip” condition on the surface of 
the plane, the flowrate out the right side of the 
system is less than the flow in the left side. The 
excess has to go out the top, causing v y > 0 there. 





=s 



=B 

3-- 



^ j 

*» 1 . 

^ 1 hJ 


6(x) 


h i I 


X 


Boundary-Layer Thickness 

The definition of boundary layer thickness is somewhat arbitrary. Although we are tempted to say 
that 8 is that value of y at which the boundary-layer solution for v x (y) equals U (potential flow), v x 
approaches U only asymptotically as y— >°o; therefore, this "definition" yields the unhelpful result that 8 is 
oo. There are several ways to assign a more meaningful finite value to 8(x). 
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One way to define the boundary-layer thickness, 8, 
is that value of y at which the velocity is 99% of the 
asymptotic value. 


U 


• Defn 1: 


v x (x,5i) = 0.99Uq(x) 


For flow over a flat plate, this convention yields: 


8 1 (.v) = 5.0, — 

iW Mu 



Another way to define 8 was suggested by Nemst who was 
concerned with boundary layers which arose in mass transfer 
problems. Nemst chose the diffusion boundary-layer thickness 
as the thickness of a hypothetical stagnant film which has the 
same diffusion resistance. 


AF 


=-iA 

v=o dy 


y = o 




Graphically, this 8 can be determined by drawing a tangent to 
the concentration profile at the surface (y=0). By analogy, 
concentration of mass is like concentration of momentum, which is just the fluid velocity. 


If we define 8 for momentum transfer in the same way, replacing concentration by v x , then: 


• Defn 2: 




dy 


A v x Uq(x)-0 

y=0 A y 5 2 -° 


For flow tangent to a flat plate, this definition yields: 


St (v) = 3.0, — 

’ Mu 


Recall that the normal component of the velocity profile was nonzero far from the plate: 


v 


lim 

y — >oo | U 


i^y) 


o.8604,PUo 
Ux 


except far downstream from the leading edge: 
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lim 

y — > oo 


v y(^y) | 

u 


lim < 0.8604 J — 


► oo 



= 0 


This means that the streamline moves away from 
the plate. A third definition for 8 is the distance 
the streamline is displaced away from the plate. 


Suppose that, far upstream, a streamline is given 
by 

a:— >-°° : y = a 



Thus the streamline is initially a distance a from x 

the a: -axis. If downstream from the leading edge 
of the plate, the streamline is a distance b from 
the A>axis, then the displacement of the streamline is given by: 


83 = b-a 

To evaluate this displacement, recall (from the definition of streamline) that the flowrate in the x- 
direction between y=0 and the streamline is the same all along the streamline. Thus: 

Q b ^ 

V = — = aU = \v x dy 

0 

where W is the width of the plate (assumed to be arbitrarily wide). Adding 83 U to both sides: 

b 

(a + d 3 )U = jv x dy + d 3 U (142) 

0 

The left-hand side of this expression can be rewritten as: 

b 

(a + 8 3 )t/ =bU =\udy (143) 

0 

b b 

J 6/ dy = J v x dy + 83 U 

(142) becomes: 0 ^ 

&3 U = \{ u ~ v x) d y 

0 
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Of course, the upper lim it b will depend on the particular streamline we have chosen. In other words, 
different streamlines are displaced differently depending on how close to the plate they were initially. 
Thus the displacement distance 83 is not unique. However, distant streamlines are all displaced by the 
same distance (since the integral converges as /?— >° o). So let's define the boundary-layer thickness as 
the displacement of external streamlines: 

00 

• Defn 3 : b 3 U = j(U -v x )dy 

0 

For flow tangent to a flat plate, this yields: 

Notice that all three of these definitions yield a boundary-layer thickness which is proportional to 
-Jvx/U although the proportionality constant varies considerably. 



We showed (see page 132) for any 2-D flow (which this is) that: 

5 ®-^== ^ 

v u 


(144) 


(145) 


where R is the radius of the cylinder. For noncircular cylinders, R is some characteristic dimension of 
the cross section (e.g. the major or minor axis of an ellipse). A semi-i nf inite flat plate is somewhat 
unusual in that it has no characteristic dimension. However, if the plate were finite with length L along 
the direction of flow, it would seem natural to choose L as the characteristic length. If one can 
reasonably assume that what happens downstream with a longer plate does not significantly effect the 
boundary layer thickness for a plate of length L (i.e. “exit effects” don’t propagate upstream). For an 
semi-infinite plate, the same result is obtained by choose x as the characteristic length; then 


8 ~ 



and 



which is consistent with (144) and with Prandtl's more general result. 


Drag on Plate 

The net force exerted by the fluid on the plate is calculated from: 
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F = \/\ n ■ I da = -\jypnda + j^n • ~gda 

Choosing n = +j on the upper surface and n = - j on the lower 
surface, the integral involving the pressure vanishes owing to 
cancellation of the contributions from the upper and lower 

surfaces (the pressure is the same, but n has opposite direction ^ 

on the two surfaces). This leaves 


2(1 


r 

dv x 


dy 

V 

y = o j 


dx 


The “2” comes from addition of the two contributions from the upper and lower surface. Evaluating the 
integral and expressing the result in dimensionless form: 

1 'x_ 

r - 2Wr _ 1-328 
D ~lpu 2 ~^ 

where Re = UxN. In defining drag coefficient this way, we have departed somewhat from the 
convention which uses the projected area along the direction of flow (i.e., the area of the shadow cast 
by the object if the light source were located very far upstream). In the case of a plate this projected 
area is zero, so we have used the area of the plate instead. 


Solution for a Symmetric Cylinder 


Let’s now return to the problem of flow around a cylinder 
at large Reynolds number. We follow the same general steps 
as we did in solving flow tangent to a flat plate: 


• Step 1 : find potential flow solution and Uq (jc) 

We accomplished this in HWK #1, Prob. 4. The velocity 
profile in the potential flow solution is 


v r =U | 


1- 


'/Cj 2 

\r J 


COS0 


v Q =-U 


1 + 


'R] 2 

K r J 


sin0 


• Step 2 : apply Prandtl’s b.l. equations (see page 134) 


r 
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dv r d\> ov, 

V — — h V V 

dx y dy c)y 2 


d2y x _rr dU 0 
- Un — ; 


dl’ 

— — H — = 0 


dx dv 


From this potential-flow solution, we can calculate the U Q appearing in Prandtl's equation: 


U 0 (x) = -vq(R,Q) = 2U sinO = 26/ sin— 

, ( X 3 , X R 

sin n =sin — 

l R ) R 


where x is the arc length, measured from the forward stagnation line, and 
0 is the polar coordinate, which is measured from the rear stagnation line. 
The two are related by 

x = R(% - 0) or 0 = 71 - 

Prandtl’s boundary-layer equations (126) become 


dv.. dv r d z v Y U z x x 

— ^ + v v — --V — f L = 4 sm— cos — 

dx y dy Oa; 2 R R R 


dx ■ dy d y 2 
dv Y dv x 


dx + dy ° 


Appropriate boundary conditions are: 
no slip: 


v x = v y = 0 at y=0 


and outside the boundary layer, we obtain potential flow: 

-» U 0 (x) as y -> oo 

The main difference between a circular cylinder and the flat plate the left-hand side of the first equation, 
which represents a pressure gradient along the direction of flow inside the boundary layer. 

Uo du SL = _}_dp #0 

dx p dx 

• Step 3 : re-write the b.l. equations in terms of the streamfunction 


In terms of the streamfunction, Prandtl’s boundary-layer equations are: 
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VyVxy -¥x¥ w ~ V V 


yy 


yyy 


. U . x x 
= 4 sin — cos — 

1 # 41 2 ^ 4 # 

dU 0 


U, 


O' 




at >’=0: = Vy = 0 

as y— > °°: \[f — > U a (x)y 

Blausius obtained the solution to this problem, which in fact works for cylinders of much more general 
shape — not just circular cylinders. 

What is required is that the potential flow solution must be an odd function of x: 

U a (x) = u | x + W3.X; 3 + M5X 5 + ... 

Then Blausius obtained a solution with the form (see S: 154-158): 

4 f(x,y) = (v/w 1 ) 1 /2 [ w 1 w/i (T| ) + 4«3.v 3 / 3 (q) + 

+ 6u$x 5 f 5 (tj ;m 1 ,m 3 ,m 5 ) + ...] 

where q = y(it\/v) l/2 

The velocity profile obtained this way for a circular cylinder is sketched below: 


This result is quite different from what was obtained with the flat plate. With the flat plate, the velocity 
profile v x as a function of y had the same shape for different v’s. Indeed the shape at different v’s were 
similar and we were able to use a “similarity transform.” Here, for flow perpendicular to a circular 
cylinder, the basic shape changes with x (or a or 0). In particular notice that the initial slope 
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@ v xM | y=0 


is positive for a=0 and becomes zero for a= 108.8°. At even larger a’s, the initial slope 


becomes negative. 


This turning around of the velocity profile coincides with a profound event called boundary -layer 
separation. Separation does not occur on the flat plate. The main difference in Prandtl’s boundary- 
layer equations which causes separation is the form of the pressure gradient. 


Boundary-Layer Separation 


To see what boundary-layer separation is and why it comes about, let's first recall the potential flow 
solution for the pressure on the surface of the cylinder. Using the velocity profile obtained in Hwk #4, 
Prob. 4, we can calculate the kinetic energy per unit volume at any point in the flow 



= U 2 !i + bl- 2 cos 2 
FE 



J_FRl 
+ 2 H r K 


4C 


Then using Bernoulli’s equation: 


we can compute the pressure profile around the cylinder. 

Consider a fluid element approaching the cylinder along 
the stagnation line shown in the sketch at right. As the 
fluid element moves toward the stagnation point A (0 = 
n and r changes from oo to R), the pressure rises to a 
maximum. Moving away from Point A along the 
surface of the cylinder ( r=R and 0 decreases from n 
to n/2), the fluid element now accelerates until reaching 
its maximum speed and minimum pressure at Point B 
and so on. 
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Hie results for the pressure changes or kinetic energy 
changes during this journey are summarized on the graph 
at right. Note that Bernoulli's equation balances an 
increase in kinetic energy with a decrease in pressure 
and vice versa. So that the sum of the two curves 
equals a constant. Inside the boundary-layer, the same 
pressure gradient applies but we also have viscous 
dissipation in addition to kinetic energy. The upshot is: 

just outside b.l . : potential flow (no imev. loss of 
energy). The kinetic energy of the fluid at B is just 
enough to overcome the pressure hill at C. Fluid 
elements arrive at C with v=0 (no kinetic energy to 
spare). 

inside b.l.: dp/dx same, but viscous dissipation 

consumes some of the kinetic energy, leaving 
insufficient energy to climb the pressure hill. 

Consequently fluid elements in the boundary layer stop 
their forward advance at some point before reaching C, which 
we will label S. Fluid elements between S and C are driven 
toward S by falling down the pressure hill. 

To conserve mass, fluid must be pushed away from the cylinder 
at S. This is known as separation of the boundary layer. 

Hie x -component of flow is toward the separation point on 
either side of S. Thus at y=0 dv x /dy> 0 for x<x s and dv x /dy< 0 
for v >,i ' s . At the separation point the derivative changes sign: 


which serves as a convenient way to locate the separation point 
in any mathematical solution to the flow problem. 

Necessary conditions for separation include: 

1. decelerating flow or Dp/Dt> 0 

2. irreversible losses of energy 
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Drag Coefficient and Behavior in the Wake of the Cylinder 

Blausius solution does not apply downstream from the separation point (i.e x>x s ). Indeed, the 
boundary layer is not thin enough to be described by Prandtl’s equations. Since the behavior on the 
downstream side significantly inflences the drag and since this behavior is not predicted by Blausius 
solution, Blausius solution does not correctly predict the drag force. Measurements of the drag 
coefficient for flow normal to a circular cylinder are summarized below. 



Several different regions are apparent in this log-log graph. Some of these correspond to major 
changes in the shape of the velocity profile. 

Re « 1. Upstream and downstream halves of streamline are Se = o.o38 

mirror images. This is what you would expect if you looked for a 
streamfunction solution with the form (see Hwk #7, Prob. 2): 

\|/(r,0) =/(r)sin0 

Although inertial terms are never negligible. The measured drag 
coefficients agree well with the prediction of Lamb. 

As we increase the Reynolds number, inertia becomes more 
important. Generally large inertia has a tendency to make 
streamlines straight. 

1 < Re < 6. Streamlines are no longer symmetric about 0=71/2. Re= 11 

Streamlines are somewhat farther from the cylinder on the 
downstream side. 

Drag coefficients are significantly smaller than predicted by Lamb. 
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6 < Re < 40. Two stationary vortexes are formed in the wake of the 
cylinder — a consequence of boundary-layer separation. 

Re > 40. Periodic vortex shedding. Large vortices, the size of the 
cylinder, are created on a periodic basis. These vortices detach and 
move downstream. Vortex shedding alternates between the top and 
bottom of the cylinder. Vortexes shed from the top have a vorticity 
(Vxv) with a sign opposite from vortexes shed from the bottom. 



Re= 140 

These vortices persist many cylinder diameters downstream from the cylinder. Most of the irreversible 
losses of energy occur in forming these vortices, whose ultimate fate is to dissipate their kinetic energy 
as heat. 

Re « 4xl0 5 . Onset of turbulence in the boundary layer. Point of separation moves further back 
toward rear stagnation point. Drag is significantly reduced — almost a discontinuous drop. 
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The Lubrication Approximation 


Consider a plate sliding on a lubricating film past a 
second stationary surface. If the distance separating the 
two plates is small compared to the dimensions of the 
plate, we can assume fully developed flow applies 
through out most of the oil film . Then: 

for h«L: v x [y) = U — - — 

p = Pq for all x,y 



L 

Q = \ v x{y) d y = \uh 
0 

^ xy ^yx 


\xU 

h 



Since the pressure is same inside the film as outside the slide, the sliding motion of two parallel surfaces 
produces no lateral component of force. 



Now suppose the slider is inclined ever so slightly 
relative to the stationary plate. We might guess, 
that if a is small enough, the velocity profile will 
not be affected. But, owing to the inclination, h is 
no longer independent of x, so our guess leads to: 

§ = ^Uh(x) = f(x) (146) 



which violates continuity. It turns out that 

continuity is preserved by a nonzero pressure gradient, dp/dx, which causes pressure-driven flow. Thus 
even the primary component of the velocity profile is affected by this slight inclination. More 
significantly, it turns out that this inclination will produce a different pressure in the film from the fluid 
outside the slider block which, tends to push the two surfaces together or apart. 

Fy* 0 
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Let’s try to estimate this force. It turns out to be no 
more difficult to obtain the result for an arbitrary gap 
profile h(x) (see figure at right), since the essential 
difficulty arises from the fact that h is not constant with 
respect to .v. 

Suppose the thickness of the gap is everywhere very 
small compared to the dimensions of the slider block. 

h(x) « L for all x 



U 


Essentially, this is a geometry with two very different 

length scales characterizing variations in the different directions x and y: we expect slow variations with 
x and rapid variations with y. We will exploit this difference using a regular perturbation in the ratio of 
the two length scales. 


We will start by nondimensionalizing the equations of motion: 


Let 




_v. 


L is an obvious choice for the characteristic value of x (since Oo c<L) and U is a obvious choice for the 
characteristic value of v x (since v x = U at y=0). Since 0<y<h(x) some characteristic value h c of the 
film thickness seems like a logical choice to scale y* The choices of characteristic values for v y and p 
are not obvious; so we will postpone a choice for now and just denote these values as v c and p c . 

We seek a solution in the form of a regular perturbation: 

where a = — 

L 

v x (x,y,a) = Uu(X,Y,a) = u[u 0 (X,Y) + aui(X,Y)+..] (147) 

v y {x,y,a) = v c v{X,Y,a) = v c [v 0 {X,Y) + a Vl {X,Y) + ..] (148) 

p{x,y,a) - Poo = Pc P( X,Y, a) = p c [P 0 (X, Y) + aP { (X , Y)+. . ] (149) 

As usual an important aspect of this form is that all derivatives of the dimensionless velocity components 
u and v with respect to the dimensionless coordinates X and Y are 0(a°): 


* For definiteness, we might select the largest value of h(x) to be h c . 
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Substituting (147)-(149) and (152) into (154): 

U 2 du aU 2 du _ p c dP U d 2 u _ _ U d 2 u 
~ 11 W + ^ v ^ ~~JlJx + v ^^ +v ( aL ) 2 ay2 

°(“°) °(“°) r_ ' o(U) 

The last term in the equation is lowest-order in the small parameter a: it’s O(or 2 ). All the other terms in 
the equation (except possibly for the pressure gradient) are 0(a°). Unless we have some other term in 

the equation of the same order, we will be forced to take ^ = w hich yields linear shear flow 

to all orders in a. We already know that this solution violates macroscopic continuity. To avoid this 
situation, we choose p c so that the pressure gradient term is also of 0(a 2 ): 
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Pc = P Lx 



| xU 
a 2 L 


With this choice, the leading term in NSE^ is 


d Uq 5/q 

~dY 2 ~ ~dX 


(155) 


(156) 


So far we have two equations [(153) and (156)] in three unknows («q, v 0 and Pq). We need another 
equation. So we turn to the secondary component of NSE: 


NSE,: 


di\ 


v Y — — + v v 
x -\ y 


c)x 


BVy 

1 dp 
£_ 

( -\2 
d Vy 

dy 

i V 

p 

v dx 2 


■ + 


-\2 A 

O Vy 


Introducing our dimensionless variables, we have: 

(. lU 


a U 2 dv (aU) dv 

u v — 

L dX a L dY 




dP a U d 2 v 


O(a) 


O(a) 


(a 2 Z,)(ocL) 


dY 


• + v- 


■ + V- 


aU 




Lr dX z (a Lf dY z 


0( a) 


O 


(a- 1 ) 


After substituting the perturbation expansions (147)-(149), (152) and (155), then collecting terms, the 
leading term is: 


a 


-3- 


dP 

-P- = 0 or Pn = const w. r. t. Y 
dY u 


(157) 


This conclusion suggests that the pressure gradient in (156) can be treated as a constant with respect to 
Y: 

d = const w.r.t. y (158) 

BY 2 dX 

We might be tempted to set this pressure gradient to zero since: 

p(0) = p(L) = p 0 

Of course the pressure gradient is zero when the two plates are parallel. But we already suspect that 
the pressure gradient is needed to avoid linear shear flow, which violates continuity. So we avoid this 
temptation and leave dP^dX nonzero. 
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From this point on, we will resort to the dimensional variables, using what we 
have learned from the leading terms of the perturbation expansion. In essence, 
we are truncating (147)-(149) after the first nonzero term. 

Integrating (158) twice: v x (x,y) = ——y 2 +c 1 (x)y+ c 2 (x) 

2(1 dx ’ 


Boundary conditions are given by the “no slip” requirement: 


v x = 0 at y=h(x ) 
v x = U at y=0 


Evaluating the two integration constants leads to: 




U 


1 - 


h 


+ 


1 dp 
2|i dx 


(y 2 - yh ) 


linear she ar flow pressure - driven flo w 


The volumetric flowrate per unit width of plate is calculated as: 


Q 

w 


h 

J v x (x,y)dy 
0 


Uh h 3 dp 

2 12p dx 


(159) 


Q can be made to be a constant at each x if dp/dx is allowed to take on non-zero values. The 
necessary values can calculated from (159): 


<2/W=const.: 


dp _ 6\xU 12|i Q 

dx h 2 h 3 W 


Integrating with respect to x from x=0 where P=Pq to any other x. 


X 

p(x) - Po = 6 pU J 
o 





(160) 


(161) 


The pressure is the same at the downstream end of the gap. Then: 


P{L ) - Po 


L 

0 = 6(l£/J 

o 




dx 

h 3 


Knowing the overall pressure drop is zero allows us to compute the flowrate: 
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where 


Q _UH 
W 2 


L r dx / L r dx 


(162) 


is an average gap width. If h(x) is a linear function whose value changes from 

h{ 0) = h\ 

to h(L) = h 2 < hi 


then 


H=2- h ' h - 


1 1 


or 


i o 

1 

[ hi h 2 ) 


hi+h 2 H 21 

which is called the “harmonic mean” of h\ and h 2 . (162) into (160) gives the pressure gradient: 
dp _ 6(lf/ ( l _H\ 


dx ip 


h 


Note that 

h 2 < H < hy 

at x=0: 

h = h 2 > H — > dp/dx > 0 

atx=x H : 

h = H — » dp/dx = 0 

at x-L : 

h = h \ < H — » dp/dx < 0 



Thus H represents not only an average gap width (with respect to flowrate) but it is also the width at the 
point where pressure is a maximum. 

Now we are in a position to evaluate the force exerted on the plate by the fluid. The x -component of 
the force is 


F x= w { V (*> 0 )dx = W j [l 


L dv , 


o dy 


,y=o 


dx = WL^ + 0(e) 
h 


where 


£ = 


h\ ~ h 2 

L 


is the angle of tilt between the two plates (| e | « 1). This result is the same (neglecting the 0(e) 
correction) as would be obtained for two parallel plates. More interesting is the v-component. For a 
linear gap: 
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J L [pO) - Po] dx = 


6| xUL z 


{hl~h 2 y 


ln- 


hn 


h\ + h 2 


Does this reduce to zero when the two surfaces are parallel (i.e. h\=h 2 )l Although it might appear that 
Fy — >oo as h | —>h 2 , it turns out that the terms inside the square brackets tend to zero faster than the 
denominator outside the brackets. In the limiting case that 


e«l : 


/q ~ h 2 ~ h. 


and 


then we obtain: 


h\-h 2 « mm(hi,h 2 ). 





L\ 


\ h J 


£ 


So we do recover zero force for parallel plates. The main difference between two parallel- and two 
nonparallel-plates is the occurance of the nonzero y-component of lift which would not occur for parallel 
plates. Notice that F y > 0 if e>0 (h\>h 2 ) and F v <0 if £<0 ih\<h 2 ). Thus either repulsion or attraction 
of the two plates is possible, depending on the direction of 
tilt relative to the direction of flow. 


Translation of a Cylinder Along a 
Plate 

The lubrication approximations developed for the slider 
block can be easily extended to other geometries. For 
example, instead of a planar slider block, suppose I try to 
drag a cylinder parallel to a plate. What will be the force 
tending to push the two surfaces apart? The same 
perturbation expansion done with the slider block applies 
here, except we have a different profile h(x) for the gap 
between the two surfaces. 

To deduce the gap profile, recall the equation of a circle 
{x-x c f + {y-y c ) 2 = R 2 

where (x c , y c ) is the location of the center of the circle and 
R is its radius. Substituting the coordinates of the center in 
our problem and y(v) = h(x), we have 



x 2 +(h-R-h 0 ) 2 = R 2 


Copyright © 2000 by Dennis C. Prieve 





06-703 


164 


Fall, 2000 


Dividing by R- 


■ + 


R- 


h - h( 
V R 


32 


= 1 


Recognizing that h-h 0 is very small compared to R, we can use the Binomial Series, truncated after the 
second term, to obtain an approximation to the second term: 


= ( 1 - e ) 2 =1 — 2 £ + d ( e 2 )- 1-2 — ^ 

R 


1 

1 

-s: 

L 

0 

1 

1 

V R ) 


^_R_^ 



l £ J 


Hie “1” cancels with the “1” on the right-hand side of the equation, leaving: 

2 


h(x) = /?() + 


2 R 


(163) 


provided that h-h 0 « R which requires that x remain small compared to R. 

dh x 

a = — = — « 1 

dx R 

where we have moved x=0 to the center of the gap, which is now symmetric about v=0. Although 
(163) is only valid for x very small, it turns out virtually all of the contribution to the force comes from 
the region where (163) is valid — provided that /i 0 is sufficiently small compared to R. In any case, let's 
assume that (163) is valid. If that bothers you, then 
replace the circular cylinder by a parabola. 


cylinder 


As with the flat slider, the pressure profile is 
determined by the need to have the volumetric flowrate 
through any ,i'=const plane be the same for all such planes. 
Eq. (159) becomes: 


Q Uh /A dp 


W 


12(1 dx 


= const w.r.t. ^ 


(164) 


If we view Q/W as an unknown integration constant, then 

we will need two boundary conditions to evaluate the two 
integration constants we will have after integrating this. 

Since the fluid held between the cylinder and the wall is in 

contact with the same reservoir at either end of the gap, we can require: 

P = Pq at X = -oo,+oo 



Copyright © 2000 by Dennis C. Prieve 




06-703 


165 


Fall, 2000 


By “infinity’’, we simply mean far from the origin. It might seem more reasonable to specify x=R, but if 
R is sufficiently large compared to h 0 , no significant error will be incurred by extending the limit to 
infinity. The counterpart to ( 1 6 1 ) is: 


Pi*) ~ Poo 



(165) 


Applying the other boundary condition at .r=+°° allow us to evaluate 

Q _UH 
W 2 


Q. The counterpart of (162) is 

(166) 


where 


is an average gap width. 
Substituting (163): 



* = 7*o 


As long as h(x) is an even function of x, then 
p(x) must be odd: 


P ~ P 0 
\lU 


h(x)=cvcn — » p(x)=odd 

Thus the pressure profile given by (165) looks 
as shown at right. The extrema correspond to: 
dp/dx = 0 

Subsisting (166) and dp!dx = 0 into (164) yields: 


dp/dx= 0: 


h = H=^h 0 



Substituting (163) and solving for v: 


dp/dx= 0: x = Rh() 

Because p is an odd function, there will be no normal force tending to separate the two surfaces: 
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/t=odd: 


F 00 

-^-= j[p(x)-p Q ]dx = 0 


W 



(*) 


dx= 271)1 f/ 



Cavitation 

However, there is an important phenomenon which we have not discussed which can cause a lift 
force to push the two surfaces apart. That phenomenon is 

cavitation - formation of gas bubbles caused by a lowering of pressure 

If the absolute pressure of the fluid drops below the vapor pressure of the liquid, we will have boiling of 
the liquid and cavitation. Because pressures generated in lubrication problems can be significant 
compared to atmospheric, cavitation is not an uncommon event. 

sources of gas bubbles : 

. vapor of liquid (if p<p vapor ) 

* air (if p<Psaturation) 

Many liquids are kept in contact with air at one atmosphere and therefore become saturated with air. If 
the pressure on the liquid is suddenly lowered, the air will be supersaturated and air bubbles will form. 

What effect will cavitation have on the pressure profile? 

Although an exact analysis would require consideration 
of two-phase flow, we can anticipate that — at the very 
least — the absolute pressure cannot drop below 
saturation. 


If any of the negative portion of the pressure profile is 
chopped off, the profile loses its anti-symmetry. A 
repulsive force pushing the surfaces apart becomes 
likely. The resulting profile might be expected to look 
something like that shown at right. 
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Squeezing Flow 


The above two problems both involve the sliding of two 
surfaces past one another. A related lubrication problem is the 
squeezing motion between two bodies. Consider the “squeezing” 
motion in a thin fi lm of liquid held between a circular disk and a 
parallel flat plate in the lim it in which h«R. In the limit e = h/R 
— » 0, a solution to the Navier-Stokes equation can be found via a 
regular pertubation expansion of the form: 

v r (r,z,£ ) = u c u{ p,C,e) = n c [n 0 (p,Q + £iq(p, C, )+...] 



w 


ft 




circular 

disk 



A 


plate 


( r,z , e) = Uv{ p, C,e) = U[v 0 (p,Q + £iq (p,C)+. • ■] 


P^iZ^-Poo = p c P(p,C,e) = p c [p 0 (p,Q + ep 1 (p,Q+...] 


where 



h_ 
R ’ 




Note by using the arbitrary u c as our characteristic radial velocity, we are delaying the choice until we 
have a chance to inspect the continuity equation: 


continuity: 


Nondimensionalizing: 


Dividing by Ullr. 


1 jK) , dv z -Q 

r dr dz 
u 1 d(pn) U d\> 

Tp dp + Jd^~ { 

£u c 1 9(pn) ^ dv 
p dp dC, 

1 o(e°) o(e°) 


( 167) 


(168) 


Boundary conditions on y _ include: 


at z= 0 or C=0: v z = 0 or v = 0 

at z=h or C=l: v z = U or v = 1 

This means that dv z /dz (and dv/dC, ) is not zero; so that the other term in the continuity equation (167) 
is not zero either. The only way the two sides of (168) can have the same order in e is for the 
coefficient to be 0(e°). This is accomplished by choosing: 

u c = E'^U 
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Next we examine the principle component of the Navier-Stokes equation. Although the flow is caused 
by motion in the z-direction, the r-component is much larger ( u c » U for 8 « 1). 


r: 


d v r 

dr 


• + v. 


dv r 

dz 


1 dp (I d 1 d(rv r ) 
Pf dr p f dr r dr 


r\ 

(I d v r 

P f dz 2 


Here we have added the subscript ” on the fluid density to avoid confusion with the dimensionless 
radial coordinate. Nondimensionalizing: 


( £ ^ du £ ~ l U 2 du _ 
R U dp + eR v a; - 

°( E ~ 2 ) °( e ' 2 ) 


dp |l£ l U d 


PfRdp p f R 2 dp 


1 d(pu) 

P dp 


+ 


(18 l U 


P /(eR) 2 3C 2 


o 


(«-) 


o 


F 5 ) 


Clearly, the first viscous term and both inertial terms are negligible compared to the second viscous 
term. Unless we have some other term in the equation of the same order, we will be forced to take 
~\2 / 

a L y 9 = 0, which (after no slip is applied) yields u - 0 for all C This would violate continuity. Thus 

/dC 

we choose p c so that the pressure derivative has the same order as the second viscous term: 


Pc l U 


or Pc = 


j xU_ 
8 3 R 


P f R P /{eR)' 

Finally, we look at the secondary component of the Navier-Stokes equation: 


z: 


dv^ dv^ 1 dp |1 1 d f dv^} |I d 2 v v 


dr 


~ + v z — 


dz 


- + ■ 


Pf dz Pf 


dr 


dr 


+ - 


Pf dz" 


Nondimensionalizing: 


V-V, 


R 


■ + v— = - 


dp 8 R dC, 


a 


(•-■) 


o\ 


dv 

/ 8 3 R dP , 

j! U Id 

f dv) 

l dp) 

LI U d 2 V 

_l_ * 


PfZR dC, 

Pf R 2 P 3p 

Pf (e«) 2 3C 2 


o( E - 4 ) 

0( E ») 


4" 2 ) 


Note that dP/d'Q is the lowest order term in this equation. Moreover, it is the only term which is 0(e~ 
4 ). Therefore when the perturbation expansions are substituted, and terms of like order are collected, 
the leading term is 
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8 ^: 


c)P< 




r = ° OT ^o = ^o(p) 


(169) 


so that the z-dependence of pressure can be neglected compared to the r-dependence. The resulting 
problem to solve (going back to dimensional quantities) is: 


Continuity: 


l3(rv r ) _ dv z _ A 

+ ^r =0 


(170) 


NSE r 




d 2 v r _ dp 
dz 2 dr 


const. w.r.t. z 


(171) 


Since the right-hand side is independent of z, we can integrate immediately to obtain the general solution: 


V- V r = + ( r ) z + c 2( r ) 


Boundary conditions are: 


at 8=0: v r =v z = 0 

at z=h : v r = 0, v z = - U 


Applying the b.c.'s we get: 

Vr =—— Z (z-h) (172) 

2 \xdr y ’ 

As before, dpldr is determined such that continuity is 
satisfied. Now, however, macroscopic continuity requires: 

h 

| v r (r, z) 2nrdz = J irty (173) 

0 flow in 

' v ' through top 

flow out through 
walls of cylinder 

(172) into (173) and requiring that the result be satisfied for 
any r yields: 


dp 

dr 



(174) 


1 

1 V V V 1 

1 

w w w > 

i 

0 

A 

3 




1 

1 

1 

r 


Side 'View 



Requiring that: 


p(R) =Po 
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we can integrate to obtain the pressure profile: 

p{r)~Po =3fit/- 


R 2 -r 2 


h 3 

which is sketched at right. Substituting (174) into (172): 

(175) 


v r (r,z ) = 3 U 


h 


r z_ \ f z ^ 2 

yh) 


\hj 



The other component of velocity can be determined by 
applying microscopic continuity. Using (151) and (175): 


c)v z _ 1 d(rv r ) 

dz r dr 


6 U 

h 


( Z^ 


\nj 


' z) 

yh) 


Integrating subject to v- = 0 at z = 0: 


v z {r,z) = -U 


U^ 2 


y h j 


f z ^ 

yhj 


To calculate the force exerted by the plate on the fluid, we use 
the unit normal pointing out of the fluid: n=k: 


O' 1 


(176) 


dF = k • T da I 

1 n — I, 

From the axisymmetry of the problem, we anticipate that there 1 

will only be a z-component of this force, which we can calculate 
by post dotting the above by k: 

dF- = k ■ T • k da = (-p+x zz )da 

In this problem, the normal component of the deviatoric stress (x zz ) vanishes. Using (176): 

x zz = p dvjdz\~h = p(-U)[6zh~ 2 - 6z 2 lr 3 \\ z=h = 0 

This leaves just a contribution from the pressure. Since p(r,h) is independent of 0, we choose da = 
(2nr)dr to be a thin annulus of radius r and thickness dr. 


R 

271 J rp(r,h)dr 
0 


37t [lUR 4 
2h 3 
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Notice that, for a fixed U. the force required becomes unbounded as h—> 0. 


Reynolds Equation 


Sliding motion and squeezing motion are quite different. Yet the lubrication approximation for each 
has a lot in common. In particular, note that in either case, pressure can be taken as as a constant along 
a normal to either surface (compare (157) and (169)). In either case, the dominant velocity component 
is tangent to the surfaces and the principle component of the NSE is approximated by a balance 
between viscous shear stress on the surface with the pressure gradient along the surface (compare (158) 
and (171)). It turns out that the lubrication approximation can be generalized to handle an arbitrary 
combination of squeezing and sliding motion in 3-D. 

Consider two bodies of arbitrary (but smooth) shape 
moving slowly through a viscous fluid in the near vicinity of 
each other. A rectangular Cartesian coordinate system is 
chosen so that the "-axis coincides with a straight line 
connecting the two surfaces at the points of minimum 
approach. The origin is located at some arbitrary point 
along this line. b i represents the distance (along the z-axis) 
from the origin to the surface of body i (i = 1,2), while z=- 
h\(x,y) and z=h 2 (x,y) describe a portion of their surfaces 
nearest the origin. Let R ix and R iy be the radii of curvature 
of body i in the x- and v-dircctions, respectively. 


v - U 2 e x + V 2 e y + W 2 e z 



For distances hj(x,y ) much less than both R lx and R jy , hj can be approximated by* 


2 2 

hj (x, y) = 8 ■ + — — + — - 
27?, 27?,- 


'IX 


iy 


The total distance between the two surfaces is 


h(x, y) = h x (x, y) + h 2 (x, y) 


8 + 


- + - 


2 R x 2 R y 


where 8 = 8 1 + 8 2 and R-. (j = x,y) can be considered to be the radii of curvature of the fi lm : 


* Actually, this assumes that the principle radii of curvature of both surfaces he either in the x- or y- 
directions. Should one of the surfaces be rotated around the z-axis by an angle 0, the function acquires 
an addition term which is proportional to xy sin0 cos0. 
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R j 


1 1 
- + ■ 


R \j R 2 j 


As in the two previous examples, scaling leads us to the following equations for the velocity profile in the 
film : 


d 2 v r dp 

I 1 f~ = -T 

dz 2 dx 

d~ v y dp 
M- ^ = 

dz 2 3y 


where the pressure profile is independent of z: 


p=p(x,y) 


Since pressure is independent of z, these equations can be easily integrated to yield the velocity profile in 
the film, which again turns out to be the sum of linear shear flow (from the sliding motion) and a 
parabolic pressure-driven flow. 


v - = U 2 ~ z(h 2 -h)~ h\h 2 ] + ^(z + hi)+U 

ax 2ll J h 


dx 2p L 


A U 
h 


(177) 


Vy = a^2p[ z2 ~ z ^ ~ h ^~ hlh2 


+ ^(z + h l )+V l 
h 


where AU = U 2 -U\ and AV = V 2 -V\ 

Still unknown is the pressure profile, which is found by requiring the velocity profile to satisfy the 
continuity equation. In particular, since pressure is independent of z, we will choose p to satisfy the 
integral of the continuity equation with respect to z 

h 2 ( x,y,t ) 

J (V- v)<rfe = 0 

(x, y, t ) 


Expanding the divergence and separating derivatives with respect to z from those with respect to x and 

y- 






( 178 ) 
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Applying Leibnitz’ rule for differentiating an integral whose lim its are functions of the differentiation 
variable: 

h 2 {x,y) h 2 (x,y ) , 

d f , x , f dv x , , , , x 5 ( -^ 1 ) 

— v^(x,y,z)Jz= - 3 — dz +v x (x,y,h 2 ) — -v x {x,y,-h l ) — 

-hi (x,y) -h\ (x,y) U 2 U\ 


SO 


hi hi 

“ 3v v . a ' 


1 

-h 


dx dz dx 


f v x (x,y, z)dz-U 2 ^-U x ^- 
J ox ox 


-h 


(179) 


Substituting the velocity profile (177) into (179) and integrating leads to: 

hi ^ 

1\ + hz (JT , TJ A (fy + ^ 2 ) 

v t< fe = — (( / 1+( / 2 ) -- 


J 


-h, 


“2 

^ | 

Differentiating — v x dz = — (Uy + U 2 ) 

ra J 2 


-/j, 


( dhy dhi ) 
+ 

d 

V dp' 

dx dx j 

dx 

12(1 dx 

K 7 


where 


h = /ij + h 2 


(179) becomes 


"2 

J 

-A, 


3v r , 1 d 

-dz = -- 


dx 12(1 civ 


^ ^ 1 dA/z 

A 6/ — — 

2 d.r 


/z 3 ^ 


v 


(180) 


where Ah = h 2 - h\ 

There is a very similar result for the integral of dv y jdy : 

h 2 


\ 

-A, 


dv 

dy 


y dz= 13 


12 (t dy 


V 3 ^ 

V 3 - v y 


I A V d -A± 

2 dy 


Finally the result for the integral of dv 7 /dz : 

h 2 h 2 

J ~/fadz= J dv z = v, (x,y,h 2 ) -v z (x,y, -h x ) = AW 


—hi 


-hi 


W 2 


W\ 


(181) 


(182) 
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Adding (180), (181) and (182), setting to zero (to satisfy (178)) and rearranging: 


dr 


h 


dp 

dx 


+ 


dy 


dp 

dy 


dA/z dA/z 

-6|lA U — — - 6|lAV —— + 12(iAlF 


dx 


dy 


(183) 


which is called Reynolds lubrication equation. The solution to this equation yields the pressure 
profile in the gap for any prescribed translation of the two bodies. Outside the gap, the pressure must 
approach the bulk pressure, which is taken to be zero 

p — > 0 as (x 2 +y 2 ) —> oo 

In the special case in which upper and lower surfaces are surfaces of revolution around the same axis, 
polar coordinates (r,0) are more convenient than (x,y) since then h \ = h\{r) and h 2 = h 2 (r). (183) can 
be written in invariant vector notation: 


V s -(fc 3 V s p) = -6|i(v 2 -v 1 ) -(n 2 - ni ) 


(184) 


where v ; - is the velocity of body i (z = lor 2) and n ; - are local normals to body i (not necessarily of unit 
length). In particular n, is defined as: 


", = V ft 

where f\{x,y,z ) = h\(x,y) + z 

and f 2 (x ,y,z) = h 2 (x ,y) - z 

That n, are local normals to body z follows from the fact that/! is a constant on surface #1 (defined as z 
= -hi) and f 2 is a constant on surface #2 (defined as z = +h 2 ). If we decompose the velocity of the two 
bodies into contributions along the z axis and in the xy plane 

V/ = v. v/ + Wjt z 


then (184) becomes 


= -6(t(vs2-'Al)-'M / T? ~h) + ]2[iAW 

sliding motion squeezing motion 


(185) 
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Example #1. Let’s reformate the sliding-flow problem for 
a plate sliding past a cylinder (see page 163). In this 
problem the equation of the lower surface (the plate) is just 

h\(x,y) = 0 

The equation of the upper surface is 

h 2 {x) = h 0 +^ 

The total gap between the two surfaces is 



2 

X 

h(x) = hi + h 2 = h 0 + - — 
2 R 


For Reynolds equation, we also need 


Hie velocity of the lower surface is 


Ah = h 2 - hi = h(x) 


U X = U Vi = 0 W, = 0 

while the upper surface is stationary: 

U 2 = 0 V 2 = 0 W 2 — 0 
The following quantities appear in Reynolds equation 


A U = -U AV = 0 AW=0 


Reynolds equation becomes 


d_ 

dx 


)i 3 — 

dx 


, TJ dh 
= 6 \xU — 

dx 


which is identical to the derivative of (164). Solving this 2 ld order ODE leads to the same pressure 
profile we determined earlier. 

Example #2. 

Now let’s reformulate the squeezing flow problem on page 167. In this problem the equation of the 
lower surface (the plate) is just 


Hie equation of the upper surface is 


hi(x,y) = 0 
h 2 (x,y) = h 
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The total gap between the two surfaces is 


h = hy + /-£ = h 


For Reynolds equation, we also need 


Ah = h 2 - hi = h 


Hie velocity of the upper surface is 


u 2 = o y 2 = o w 2 = -u 


while the lower surface is stationary: 


Ui = 0 Vi = 0 W y = 0 

The following quantities appear in Reynolds equation becomes 


AU = 0 AV = 0 AW=-U 


(185) becomes 


V s ■ h 3 V s p) = -12 \xU 


(186) 


Because the upper surface is a circular disk and the gap is uniform, we expect squeezing flow to be 
axisymmetric in cylindrical coordinates. In other words, we expect that p = p(r) (i.e. no 0- 
dependence). In cylindrical (r,0,z) or polar coordinates (r,0), the gradient is* 


while the divergence is 


(186) becomes 


V p = — e 

v sr c r 


dp 

dr 




1 d ( dp ' 
dr 


1 d ^ rh 3 dp ^ 


r dr 


dr 


-12 [lU 


Multiplying through by r and integrating: 


rh 3 ^ = -6 pUr 2 + c 
dr 


* see http://www.andrcw.cmu.cdu/coursc/06-703/V ops cyl.pdf 
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Dividing by r/z 3 : 


dp r tt r 

— - -6[iU — + 


dr 


Ir 


c 

rh 3 


When this is integrated a second time, the second term will lead to a logarithmic singularity at r= 0. To 
keep pressure finite, we choose c=0 and then the above equation is identical to (174). Solving this 2 nd 
order ODE leads to the same pressure profile we determined earlier. 


Example #3. Sliding of a plate past a sphere 

n this problem the equation of the lower surface (the plate) 
is just 

h\(x,y) = 0 

The equation of the upper surface is 

r 2 

h 2 {r) = h 0 + — 

The total gap between the two surfaces is 



h(r) = h +h 9 = hn + - 


2 R 


For Reynolds equation, we also need 


Ah = /?2 - h\ = h(r ) 


The velocity of the sphere is purely along the v-axis 


v 2 = y s2 = = e r UcosQ - e e U sin0 

u 2 = u V 2 = 0 W 2 = 0 

while the lower surface is stationary: V] = v ?1 = 0 

The following quantities appear in (1 85): 

dh 

VAh-l\) = V,h = —er 

dr 


( v ^2 - V ^t) v ^ { h 2 - h \)={ u cos0 e r ~U sin 6 e e ) ■ 


dh 

— ( 

dr 


tt dh a 

U — cos0 
dr 
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(185) becomes 


ld_ 

r dr 



1 d 

' h 3 dp ' 

{ dr) 

r 30 

r dO 

K 7 


TT dh n 

= -6(1 U — cos 0 
dr 


The solution for sphere of radius R. sliding along a plate at speed U in the x -direction, is: 


p{x,y) 


6\lUx 

5[h(x,y)f 


This produces no force in the z-dircction (pressure profile is antisymmetric) but of course a force must 
be applied to the sphere to get it to move: 

F x = y 7tp UR In + o(8° ) as 8-> 0 

F y = F z = 0 
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Turbulence 

General Nature of Turbulence 

In all the problems we have analyzed to date, the fluid elements travel along smooth predictable 
trajectories. This state of affairs is called: 

laminar flow - fluid elements travel along smooth deterministic trajectories 

These trajectories are straight parallel lines for simple pipe flows. However, this is not the only solution 
to the equations of motion. Consider the following experiment 

Reynolds Experiment (1882) - inject a thin stream of dye into a fully developed flow in a pipe; 
observe the dye downstream, (see S:37) 

Flow pig 32. Water, upper: velocity 1 1 cm s, Re - 1.5 x 10 J , middle: velocity 17cnvs, Re 2.34 p 
Direction X 10 J , lower: velocity 54 cm s. Rc - 7.5 X 10', pipe ID 14 mm, dye injection method. 

1500 
2340 

7500 



• for laminar flow: dye stream appears as a straight colored thread 

As the total flow rate of fluid in the pipe is increased, a sudden change in the appearance of the dye 
stream occurs. The thread of dye becomes more radially mixed with the fluid and, far enough 
downstream, its outline becomes blurred. 

• for turbulent flow: irregular radial fluctuations of dye thread 

Using a pipe with a sharp-edge entrance, Reynolds determined the critical flow rate for a large number 
of fluids and pipe sizes. He found in all cases, the transition occurred at a critical value of a 
dimensionless group: 


P(v 2 )D 



2300 ±200 


Copyright © 2000 by Dennis C. Prieve 


06-703 


180 


Fall, 2000 


where (v z ) is the cross-sectional average velocity (= volumetric flowrate / pipe area). Today, we 
know this dimensionless group as the Reynolds number. 

origin of turbulence - instability of laminar-flow solution to N-S eqns 

instability - small perturbations (caused by vibration, etc.) grow rather than decay with time. 

That the laminar-flow solution is metastable for Re>2100 can be seen from Reynolds experiment 
performed with a pipe in which disturbances are minimized: 

• reduce vibration 

• fluid enters pipe smoothly 

• smooth pipe wall 

Under such conditions, laminar flow can be seem to persist up to Re = 10^. However, just adding 
some vibrations (disturbance) can reduce the critical Re to 2100. The onset of turbulence causes a 
number of profound changes in the nature of the flow: 

• dye thread breaks up — streamlines appear contorted and random 

• sudden increase in Ap/L 

• local v z fluctuates wildly with time 

• similar fluctuations occur in v r and v e 

As a consequence of these changes, no simplification of the N-S equation is possible: v r , v 0 , v z and p 
all depend on r, 0, z and t. 


Turbulent Flow in Pipes 

Velocity profiles are often measured with a pitot tube , which is a device with a very slow response 
time. As a consequence of this slow response time, the rapid fluctuations with time tend to average out. 
In the descriptions which follow, we will partition the instantaneous velocity v into a time-averaged value 
v (denoted by the overbar) and a fluctuation v' (denoted by the prime): 

v = v + v' 

Cross-sectional area averages will be denoted by enclosing the symbol inside carets: 
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R 

J v z (r)2nrdr 

0 

R 

J 2 nr dr 
0 


Q 

nR 2 


In laminar flow, the velocity profile for fully developed flow is parabolic in shape with a maximum 
velocity occurring at the pipe center that is twice the 


cross-sectional mean velocity: 

In turbulent flow, the time-averaged velocity profile 
has a flatter shape. Indeed as the Reynolds number 
increases the shape changes such that the profile v z 
becomes even flatter. The profile can be fit to the ( V / ) 
following empirical equation: 

^O) =^,max 




where the value of the parameter n depends on Re: 


Re = 

4.0xl0 3 

2.3xl0 4 

l.lxlO 5 

l.lxlO 6 

2.0xl0 6 

3.2xl0 6 

n = 

6.0 

6.6 

7.0 

8.8 

10 

10 

v m.a.xJ <v 7 > 

1.26 

1.24 

1.22 

1.18 

1.16 

1.16 


Hie reduction in the ratio of maximum to 
average velocity reflects the flattening of the 
profile as n becomes larger. Of course, this 
equation gives a “kink” in the profile at r=0 
and predicts infinite slope at r=R. so it 
shouldn’t be applied too close to either 
boundary although it gives a reasonable fit 
otherwise. 

How big are the fluctuations relative to the 
max im um velocity? Instantaneous speeds can 
be obtained for air flows using a hot-wire 
anemometer . This is simply a very thin wire 
which is electrically heated above ambient by 
passing a current through it. As a result of 
electrical heating (fiR) the temperature of the 
wire will depend on the heat transfer 
coefficient, which in turn depends on the 



0.0 0.2 0.4 0.6 0.8 1.0 


R-r 

R 
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velocity of flow over the wire: 


v z T 


-> 


h T 


_> T w ire T a ir ^ 


It’s easy to determine the temperature of the wire from its electrical resistance, which generally increases 
with temperature. The reason for making the wire very thin is to decrease its thermal inertia. Very thin 
wires can respond rapidly to the rapid turbulent fluctuations in v z . 


Anyway, the in stantaneous speed can be measured. Then the 
time-averaged speed and the fluctuations can be calculated. The 
root-mean-square fluctuations depend on radial position, as 
shown at right. Typically the axial fluctuations are less than 10% 
of the maximum velocity whereas the radial fluctuations are 
perhaps half of the axial. 

Note that the fluctuations tend to vanish at the wall. This is a 
result of no-slip (applies even in turbulent flow) which requires 
that the instantaneous velocity must vanish at the wall for all time, 
which implies that the time average and the instantaneous 
fluctuations must vanish. 



Time-Smoothing 

As we will see shortly, these fluctuations profoundly increase transport rates for heat, mass, and 
momentum. However, in some applications, we would be content to predict the time-averaged velocity 
profile. So let’s try to time-average the Navier-Stokes equations with the hope that the fluctuations will 
average to zero. 

First, we need to define what we mean by a time-averaged quantity. Suppose we have some property 
lik e velocity or pressure which fluctuates with time: 


s = s(t) 

We can average over some time interval of half width At: 


s(t) 


1 


j; 


2 At J t-At 


s(t')dt' 


We allow that the time-averaged quantity might still 
depend on time, but we have averaged out the rapid 
fluctuations due to turbulence. 



Now let’s define another quantity called the fluctuation 
about the mean: 
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s'(t ) = s(t ) - s(t ) 


Time-Smoothing of Continuity Equation 

The simple functional form of our experimentally measured velocity profile — V z (r) — is exactly the 
same as for laminar flow. This suggests, that if we are willing to settle for the time-averaged velocity 
profile, then I might be able to get the result from the NSE. Let’s try to time-smooth the equation of 
motion and see what happens. We will start with the equation of continuity for an incompressible flow: 


V • v = 0 


Integrating the continuity equation for an incompressible fluid and dividing by 2 At: 


1 rt + Af , , 1 rt + At 

V ■ vdt = 

2At J t-At 2/At It- At 


0 df = 0 


Thus the right-hand-side of the equation remains zero. Let’s take a closer look at the left-hand side. 
Interchanging the order of differentiation and integration: 


rt + At \ 1 rt + At I 

V-\dt' = V-\ ydt'\ = V - v 

Jt— At I 2 At t—At ! 


2At J t-At 


Substituting this result for the left-hand side of the continuity equation, leaves: 

V- v = 0 


Hius the form of the continuity equation has not changed as a result of time-smoothing. 


Time-Smoothing of the Navier-Stokes Equation 

Encouraged by this simplification, we try to time-smooth the Navier-Stokes equation: 

p^- + pv ■ Vv = -V/?+ |iV 2 v + pg 
at 

After integrating both sides with respect to time and dividing by 2At, we can break the integral of the 
sum into the sum of the integrals. Most of the terms transform in much the same way as the left-hand 
side of the continuity equation. The result is 

dv 9 

p — + pv- Vv = -Vp + (iV v + pg 
at 

With a little additional massaging (see Whitaker), the remaining term can be expressed as 
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pv ■ Vv = pv ■ Vv + V ■ (pv V) 

If the second term on the right-hand side were zero, then NSE after time-smoothing would have exactly 
the same form as before time-smoothing. Unfortunately, this term is not zero. Although the average of 
the fluctuations is zero, the average of the square of the fluctuations is not zero. So this second term 
cannot be dropped. Thus the time-smoothed Navier-Stokes equation becomes: 

p-^ + pv- Vv = -Vp +(iV 2 v + pg + V-x (r) 
dt 

where x^ = -pv / v / 

has units of stress or pressure and is called the Reynold’s stress. Sometimes it is also called the 
turbulent stress to emphasize that arises from the turbulent nature of the flow. The existence of this 
new term is why even the time-averaged velocity profile in side the pipe is different from that during 
laminar flow. Of course, our empirical equation for the V z (r ) is also different from that for laminar flow. 

Although we don’t yet know how to evaluate this Reynolds stress, we can add it to the viscous stress 
and obtain a differential equation for their sum which we can solve for the simple case of pipe flow. 
Here’s how we do it. First, recall that for incompressible Newtonian fluid, the stress is related to the 
rate of strain by Newton’s law of viscosity. Time smoothing this constitutive equation yields: 

x = (I 

_ 9 

Taking the divergence: V ■ x = jiV " v 

If we now make this substitution, NSE becomes 

p-^ + pv-Vv = -Vp + V-x + pg + V-x (r) 
dt = (165) 

= -Vp + pg + V -x (r) 

where x^- 1 =x + x (?) 

is the total stress, i.e., the sum of the time-averaged viscous stress and the Reynolds stress. Thus we 
see that the Reynolds stress appears in the equations of motion in the same manner as the viscous stress. 
Indeed the sum of the two contributions plays the same role in turbulent flows that the viscous friction 
played in laminar flow . 


Vv + (Vv) r 
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Analysis of Turbulent Flow in Pipes 


We can make the same assumptions (i.e. the same guess) about the functional form of the time- 
averaged velocity and pressure profile in turbulent flow that we made for laminar flow: we will assume 
that the time-averaged velocity profile is axisymmetric (v e =0, 3/30=0) and fully developed (3/3z=0). 


v z =v,(r) v r = v e = 0 
P = P(z) 


Then the z-component of (165) yields: 


3 P 

jk' 

f(z) 


1 3 I 


0 = — — 1- ——I rx 

r or v 


(T) 

rz 


)- 


S{ r ) 


i 4 r) 

r 

0 




dz 

IT 


where P is the time-averaged dynamic pressure. This form of this equation was obtained using the 
tables in BSL (top half of p85, eqn C), after replacing the instantaneous quantities by their time 
averages, except that the instantaneous viscous stresses x has been replaced by (minus) the total stress 
x(T> Expecting the time-averaged flow to be axisymmetric (3/30 = 0) and fully developed (3/3 z = 0, 
except for pressure), the last two terms in this equation can be dropped and the second term is a 
function of r only. This leaves us with the same equation we had for laminar flow: a function of r only 
equal to a function of z only. The only way these two terms can sum to zero for all r and z is if both 
equal a spatial constant: 


dz r dr ' rz ' J 




AP 


<0 


This implies that pressure P varies linearly with z. Solving for the total stress xj.P by integrating: 


_ 

_ 


j_A P 
2 ~L 


-r H — — x rz + x r7 


( t ) 


(166) 


The integration constant c was chosen to be zero to avoid having the stress unbounded at r=0. Now 
this is the total stress: the sum of the Reynolds stress 


r(0 __ 


“rz 


P v 'rV' z 
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and a viscous contribution from time-smoothing Newton’s law of viscosity: 

dv 7 

X rz = F— ^ 
dr 

The latter can be determined by differentiating the time-averaged 
velocity profile. If we subtract this from the total we can determine 
x\’J — one of the components of the Reynolds stress tensor. The 
result is shown in the figure at right. Notice that the turbulent stress 
tends to vanish near the wall. This can be explained by noting that at 
the wall, “no slip” between the fluid and the stationary wall requires 
that the instanteous velocity, as well as its time average, must be zero: 

^=^=0=>^=0=>xg ) =-p^ = 0 

In terms of the relative importance of these two contributions to the total, one can define three regions: 

1. turbulent core : tW»t. This covers most of the cross section of the pipe. 

2. laminar sublayer : ri'kxT. Very near the wall, the fluctuations must vanish (along with the 
Reynolds stress) but the viscous stress are largest. 

3. transition zone: Neither completely dominates the other. 

When applied to the situation of fully developed pipe flows, continuity is automatically satisfied and the 
time-smoothed Navier-Stokes equations yields only one equation in 2 unknowns: 

2 unknowns: v z (r) and pv'v' 



Clearly another relationship is needed to complete the model. This missing relationship is the 
constitutive equation relating the Reynolds stress to the time-smoothed velocity profile. One might 
be tempted to define a quantity like the viscosity to relate stress to the time-averaged velocity. 




dv z 

dr 


But if you define the “turbulent viscosity” this way, its value turns out to depend strongly on position. 


P— = {~i oo 

p [0 


near pipe centerline 
at pipe wall 


So un lik e the usual viscosity, is not a material property (since it depends on position rather than 
just the material). 
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Prandtl’s Mixing Length Theory 

The first successful constitutive equation for turbulence was posed by Prandtl in 1925. Prandtl imagined 
that the fluctuations in instantaneous fluid velocity at some fixed point were caused by eddies of fluid 
which migrate across the flow from regions having higher or lower time-averaged velocity. 

eddy - a packet of fluid (much larger than a fluid element) which can undergo random migration across 
streamlines of the time-smoothed velocity field. 

These eddies have a longitudinal velocity which corresponds to the time-average velocity at their 
previous location. 

As this eddy moves across the streamlines, it 
gradually exchanges momentum with the surrounding - v 
fluid which is moving at a different longitudinal 
velocity. But this exchange does not occur 
instantaneously. The eddy retains its original velocity 
for a brief period of time. We might call this the 
mixing time. During this time, the eddy migrates 
laterally a distance / called the mixing length: 

mixing length (/) - characteristic distance an eddy 
migrates normal to the main flow before mixing 

Although momentum exchange between eddies occurs continuously in actual turbulent flow, Prandtl 
imagined that a migrating eddy keeps all of its original velocity until it migrated a distance / and then 
suddenly it exchanges it. This is like molecules of a gas retaining its momentum until it collides with 
another gas molecule, which causes a sudden exchange of momentum. Indeed, you might find it helpful 
to think of the mixing length as being the analogue of mean-free-path in the kinetic theory of gases. 
Recall that: 

mean-free path - average distance a gas molecule 
travels before colliding with another gas molecule. 

Now suppose we are monitoring the instantaneous 
velocity at a distance y from the wall when an eddy drifts 
into our location from y+l. Because this migrating eddy 
has a higher velocity than the average fluid at y, we will 
observe an positive fluctuation when the eddy arrives. 

To estimate the magnitude of the fluctuation, we can 
expand the time-smoothed velocity profile in Taylor 
series about y=y: 



eddy migrating 
randomly 
across flow 


v x (y) 
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— 2 — 

_ . .. _ . . dv x , 1 d v x ,2 

v x (y + 0 - v x(y) + — t~ t + w i +■■■ 

dy v 2 dx- 

y y 

Assuming that / is sufficiently small that we can truncate this series without introducing significant error: 

(<\b m e =W * ( - y + l) - V * (y ' , = 1 ^- 

where the subscript “above” is appended to remind us that is the fluctuation resulting from an eddy 
migrating from above. At some later time, another eddy might migrate to our location from below, 
producing a negative fluctuation in velocity: 

« \elow = ^ y - « - 


Of course the average fluctuation is zero: v x = 0 , but the average of the squares is not: 



/ ,\2 

+ 

/ ,\2 ! 

Vx ) 

Vx ) 

above 


below J 



(167) 


Now let’s turn our attention to vC . This is related to how fast the eddies migrate, and the sign depends 
on whether they are migrating upward or downward. 

If the eddy migrates from above, it represents a 
negative y-fluctuation (it is moving in the -y direction), y 
Such an eddy will have a greater x -velocity than the 
fluid receiving it, consequently generating a positive x- 
fluctuation: 

V y ’ <0 — > V x ’ >0 — » V x ’ Vy ’ <0 

On the other hand, if the eddy migrates from below, it 
represents a positive y-fluctuation but has less x- 
velocity than the fluid receiving it, generating a negative x -fluctuation: 

Wy’>0— >V^-’<0 — > V X ’ Vy ’ < 0 

Finally, if there is no vertical migration of eddies, there is no reason for the x -velocity to fluctuate: 

v y ’ = 0 v x ’ = 0 

These three statements suggest that the y-fluctuations are proportional to the ^-fluctuations, with a 
negative proportionality constant: 
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Vy ~ av^ 


where a>0. Alternatively, we can write: 

v x 'vy' =-a(v x ') 2 
Time averaging and then substituting (167): 


v' x v'y = -a (v' x r = -a V 


v d y > 


Absorbing the unknown a into the (still unknown) mixing length parameter: 

f M 77 \ 2 


= -PVxVy = 

Comparing this result with Newton’s law of viscosity: 


dv x 

v d y y 


^ xy ~ M-' 


<&x 

dy 


we could conclude that an apparent turbulent viscosity is given by: 


|l (r) = p l 2 


dy 


(168) 


Of course, this viscosity is not a hue fluid property, because it depends strongly on the velocity profile. 

For this theory to be useful, we need a value for the “mixing length” l. There are two properties of / 
which we can easily deduce. First of all, / was defined as the distance normal to the wall which the 
eddy travels before becoming mixed with local fluid. Clearly, this mixing must occur before the eddy 
“bumps” into the wall, so: 

Property #1: l<y 

where y is the distance from the wall. Secondly, we know from no-slip that the fluctuations all vanish at 
the wall. Consequently, the Reynolds stress must vanish at the wall. Since the velocity gradient does 
not vanish, we must require that the mixing length vanish at the wall: 

Property #2: 1=0 at y=0 

If it’s not a constant, the next simplest functional relationship between / and y which satisfies both these 
properties is: 
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1 = ay 


where a is some constant and 0<a< 1 . 


Prandtl’s “Universal” Velocity Profile 


The velocity profile in turbulent flow is essentially flat, except near the wall where the velocity gradients 
are steep. Focussing attention on this region near the flow, Prandtl tried to deduce the form for the 
velocity profile in turbulent flow. Recall from (166) that in pipe flow, the total stress varies linearly from 
0 at the center line to a maximum value at the wall: 



1 A P 
2~L 


T 0 — <0 
U R 


(170) 


where we have defined xq = -( 1 /2)(R/L)AP > 0 


which represents the stress on the wall. In the “turbulent core”, the Reynolds stress dominates the 
“laminar” stress; then substituting (168) through (170): 


-KO ~ r<n 
L xy L rz 



(171) 


(172) 


where C is the integration constant, v* is called the friction velocity and where we have introduced 
dimensionless variables: 



v 


, V * 

and R + = — R 

v 


Copyright © 2000 by Dennis C. Prieve 



06-703 


191 


Fall, 2000 


Near the wall (i.e. for y«R or v + «/C~), we can simplify (172): 


for y + «R + : 


1 - 


R 


l-^— + 0(y +2 ) 
2 R + 


anh -'Jl-yl = Ll n ^l + o(y*) 

y 


r+ 2 “ -+ 


Dropping the higher-order terms: 


( +\ 2-ln(4i? + ) 1, + 

[y ) = + C — lny 


(173) 


(174) 


where c is a collection of constants. 

This result can be derived more easily by starting over with a simplified (171): 
for y«R: 


dv x . dv x 1 dy 1 dy + 

ay — - = v * or — - = — = — 

dy ay a y + 


d\- 


which integrates to 


1 


v + = — lny + + c 


(175) 


where c is some integration constant. When 
plotted on semi-log coordinates, experimental 
velocity profiles do indeed show a linear 
region which extends over a couple of 
decades of y + values: 

Moreover, the slope and intercept of this 
straight line don't seem to depend on the 
Reynolds number. Indeed, the slope and 
intercept also don't seem to depend on the 
shape of the conduit. Rectangular conduits 
yields the same velocity profile on these 
coordinates. This is called Prandtl’s 
Universal Velocity Profile : 



y + >26: 


v + =2.51ny + +5.5 


(176) 
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which applies for y + >26 (the turbulent core). This coefficient of lnv + corresponds to <a=0.4. so (169) 
becomes: 


/ = 0.4 y 


Recall that we reasoned that a had to be in the range of 0 to 1 to be physically realistic. 

In the laminar sublayer , Reynolds stress can be totally neglected, leaving just viscous stress. This 
close the wall, the total stress is practically a constant equal to the wall shear stress t 0 : 


y « R: 




xy 

dv 


fCO 

l rz 


dy 


- = *o 


Then we can integrate the above ODE for v x , subject to v x = 0 at y=0 (i.e. no slip): 


x 0 

Vj = — y 

E 


Dividing both sides by v * we can make the result dimensionless: 
y + <5: 


v + = y + 


(177) 


which applies for 0<y + <5 (the laminar sublayer). Of course (176) also does not apply near the 
center of the pipe, since the y + ~ R + there, whereas (176) was derived by assuming that _y+ « R + (see 
(173)). 


Prandtl’s Universal Law of Friction 

Let’s try to figure deduce the analog of Poisueille's Formula (see page 86) for turbulent flow. 
Poisueille's Formula is the relationship between volumetric flowrate through the pipe and pressure drop. 
Volumetric flowrate Q is calculated by integrating the axial component of fluid velocity of the cross 
section of the pipe: 


/- \ Q 2 rR — 

{ v z) = V = ^ L rv z( r )dr 

x ' HR 2 R 2 Jo 

Now we are going to use (176) for the velocity profile, although we assumed in (175) that y«R (where 

y = R-r). 
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Although the derivation of (176) assumed that 
we are very close to the wall (y+«/?+, see 
(173)), (176) works remarkly well right out 
to the centerline y + =R + . The plot at right 
shows the velocity profiles (with different wall 
roughness on the walls) compared with 
predictions based on (176). The ordinate is 


v z (R)-v z (r) _ v+ , R+ 


)-v + (y + ) = [2.5h 
2.5(ln/? + -In y + j = 2.51n 


2.51n — = 2.51n ■ R 


R-r 


Note that (176) predicts an infinite velocity- 
difference at y=0, whereas the actual velocity 
must be finite. Of course, (176) does not 
apply right up to the wall because very near 
the wall the Reynolds stresses are not 
dominant. 



Substituting (176) and integrating: 




(v*R^ 


2.5 In 

+ 1.75 


l v J 



(178) 


Now the friction velocity can be related to the friction factor, whose usual definition can be expressed 
in terms of the variables in this analysis: 


Thus 


/- 



f 


\ 2 


= 2 


v 



J 



Likewise, the usual definition of Reynolds number yields 
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Thus 


Re = 



v 


v* R 

v 


(vz)R 

-i X 




Re 77 

2 a/2 


Relating v* to /and (v z ) to Re, (178) can be written as: 


or 


77 


1.77 ln(Re 77) -0.60 




4.07 log 10 (Re 77) -0.60 


which fits experimental data remarkable well. A slightly better fit can be obtained by adjusting the 
coefficients: 


77 


4.0 log 10 (Re 77) -0.40 


(179) 


which is called Prandtl’s ( universal ) law of friction. It applies virtually over the entire range of 
Reynolds numbers normally encountered for turbulent pipe flow: 2100 < Re < 5 x KA 



Copyright © 2000 by Dennis C. Prieve 



